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FOREWORD 


NASTRAN (NASA STRUCTURAL ANALYSIS) is a large, comprehensive, 
nonproprietary, general purpose finite element computer code for 
structural analysis which was developed under NASA sponsorship 
and became available to the public in late 1970. It can be obtained 
through COSMIC (Computer Software Management and Information Center) , 
Athens, Georgia, and is widely used by NASA, other government 
agencies, and industry. 

NASA currently provides continuing maintenance of NASTRAN 
through COSMIC. Because of the widespread interest in NASTRAN, 
and finite element methods in general, the Eighth NASTRAN Users' 
Colloquium was organized and held at the Goddard Space Flight 
Center, October 30-31, 1979. (Papers from previous colloquia held 
in 1971, 1972, 1973, 1975, 1976, 1977, and 1978 are published in 
NASA Technical Memorandums X-2378, X-2637, X-2893, X-3278, X-3428, 
and NASA Conference Publications 2018 and 2062.) The Eighth 
Colloquium provides some comprehensive general papers on the 
application of finite element methods in engineering, comparisons 
with other approaches, unique applications, pre- and post-processing 
or auxiliary programs, and new methods of analysis with NASTRAN. 

Individuals actively engaged in the use of finite elements 
or NASTRAN were invited to prepare papers for presentation at 
the colloquium. These papers are included in this volume. No 
editorial review was provided by NASA or COSMIC, however, detailed 
instructions were provided each author to achieve reasonably con- 
sistent paper format and content. The opinions and data presented 
are the sole responsibility of the authors and their respective 
organizations . 

Cochairmen : 

Robert L. Brugh 

NASTRAN Project Manager 

COSMIC 

University of Georgia 

Athens, Georgia 

and 

Reginald S. Mitchell 

Goddard Space Flight Center 

Greenbelt, Maryland 
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ON THE CONGRUENT FEATURE IN NASTRAN 


P. R. Pamidi 

Computer Sciences Corporation 
Silver Spring, Maryland 


ABSTRACT 


The congruent feature, which is a capability in NASTRAN that can contribute to 
significant increases in computational efficiencies, is discussed in this paper. The 
usage of the capability and the software design characteristics affecting it are explained. 
The factors affecting the efficiency of the feature are pointed out. The details per- 
taining to the software design of the congruent feature are presented; in particular, the 
congruent element table is described. Several examples employing the congruent fea- 
ture are considered and comparisons of EMG (Element Matrix Generator) module CPU 
times with and without this feature are presented. The results of the paper clearly 
demonstrate the role of the congruent feature in increasing computational efficiencies 
and its applicability to large-size problems. 


INTRODUCTION 


An important step in any NASTRAN problem is the generation of element matrices 
(stiffness, mass, and damping matrices as required) in the EMG module. In many cases, 
this step can represent a significant portion of the total problem activity. Because of 
the differences in algorithms and procedures, the cost of generating the element matrices 
for an element depends on the element type, its configuration and its properties. How- 
ever, this cost is associated primarily with CPU activity and is not significantly affected 
by core size or I/O transfers (Reference 1). 

Normally, the element matrices are generated in the EMG module once for each 
element in the model. However, when two or more elements in the model have the same 
element matrices, there is no reason why the same matrices should be computed sep- 
arately for each such identical element. By declaring such elements as congruent, it is 
possible to cause their element matrices to be computed only once for all elements in the 
congruent set instead of their being computed repeatedly for each of the individual 
elements in the set. This results, in general, in a saving of CPU time in the EMG 
module. In many cases, judicious formulation of the problem to facilitate the use of the 
congruent feature can result in substantial savings in the computational effort. In some 
problems, over 99 percent reductions in EMG module CPU times have been obtained. 


The congruent feature is not yet adequately publicized in the NASTRAN docu- 
mentation. Currently, it is only referenced in the NASTRAN User's Manual (Refer- 
ence 2) with a one-page description of the CNGRNT bulk data card. It is hoped that the 
discussion of this feature herein will lead to more widespread use of this capability in 
large-size problems thus resulting in significant increases in computational efficiencies. 


CONGRUENT FEATURE USAGE 


The congruent feature is specified in NASTRAN by means of one or more CNGRNT 
cards in the Bulk Data Deck (Reference 2). Any number of such cards may be employed. 

The CNGRNT bulk data card is an open-ended card and requires the specification 
of a primary element identification number and one or more secondary element identif- 
ication numbers. (A description of the CNGRNT card is given in the Appendix.) The 
terms primary and secondary as used with regard to congruent data are purely relative 
and have no real significance. Generally, the primary element is the lowest numbered 
element in the congruent set, but this need not be so. The element matrices are 
computed in the EMG module only for the lowest numbered element in a congruent set 
(even though this element may not be the primary element). The element matrices for 
the rest of the elements in the congruent set are then derived from these computed 
matrices. 


SOFTWARE DESIGN CHARACTERISTICS AFFECTING 
CONGRUENT FEATURE USAGE 


When using CNGRNT cards, the user should be aware of the following important 
characteristics of the congruent capability software design in NASTRAN. 

• User Responsibility for Congruency Specification 

The elements declared as congruent must have characteristics (such as their 
orientation and geometry) that cause their element matrices in the global coordinate 
system to be truly identical. The program cannot test the validity of this structural 
specification. It is, therefore, the user's responsibility to ensure that element 
congruence specifications are valid. Improper congruence specifications will result 
in an improper structure definition and will in turn lead to erroneous results. It should 
be emphasized that the proper use of the congruent feature will not cause the answers to 
be any different from those obtained without the use of the feature, but will result in a 
saving of CPU time in the EMG module. 
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• Flexibility in Specifying Congruencies 

Clearly, congruency by its very definition can apply only to elements of the same 
type. Thus, for instance, a bar element can be congruent only to another bar element 
and not to a plate element. However, because of the effective manner in which the 
congruent feature has been incorporated into NASTRAN (as will be evident from the 
discussion in a later section), elements of different types can be specified on the same 
logical CNGRNT card without in any way making the different element types congruent. 
Thus, on the same logical CNGRNT card, several bar elements can be declared as 
belonging to a congruent set and several plate elements can be specified as belonging 
to a separate congruent set. However, the user should ensure that such specifications 
do not lead to erroneous declarations when elements of different types have the same 
identification numbers. 

• Provision of "Phantom" Element Identification Numbers 

As a corollary to the above, it may be noted that the element identification 
numbers (primary or secondary) specified on a CNGRNT card need not all exist in a 
model. This facilitates the use of the THRU option on the card more often than is pos- 
sible in many other similar cases. 

• Primary Element Specification 

The same element can appear as the primary ID on more than one CNGRNT card, 
but an element listed as a primary ID on one CNGRNT card cannot be listed as a 
secondary ID on another CNGRNT card. However, if a primary ID is also listed as a 
secondary ID on the same card, then such secondary IDs are ignored. 

• Secondary Element Specification 

The same secondary ID cannot be listed as congruent to two or more different 
primary IDs. 

• Redundant Specifications 

Redundant specifications on CNGRNT cards are ignored. 


FACTORS AFFECTING CONGRUENT FEATURE EFFICIENCY 


As indicated earlier, the use of the congruent feature results in increased 
computational efficiency. The degree of efficiency obtained depends on the following 
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factors some of which can be influenced by the user input specifications: 

• Number of Congruent Elements 

Clearly, the larger the number of elements in a congruent set and the larger the 
number of sets, the higher the savings in CPU time. 

• Type of Elements Specified as Congruent 

Larger savings in CPU time are obtained for certain element types than for other 
element types. Thus, for instance, declaring two IHEX3 elements as congruent will 
result in more savings than declaring two IHEX1 elements as congruent. 

• Type of Analysis 

For a specified congruent set, larger savings are obtained in dynamic analysis 
than in static analysis since, in the former, mass and/or damping matrices need to be 
computed in addition to stiffness matrices. 

• Numbering of Grid Points of the Congruent Elements 

Processing is slightly more efficient if the relative order of the numbering of the 
grid points of the congruent elements is the same. Thus, for instance, two congruent 
quadrilateral plate elements are processed more efficiently if their grid points are 
numbered 1-7-4-6 and 12-23-16-20, respectively, than if they were numbered 1-7-4-6 
and 11-14-17-15, respectively. In the former case, the grid point numbers of the two 
congruent elements increase or decrease in the same order as we go around the elements. 
In the latter case, the grid point numbers of the two congruent elements increase or 
decrease in different orders as we go around the elements. 


SOFTWARE DESIGN OF THE CONGRUENT FEATURE 


The preliminary checking of the validity of the data on the CNGRNT bulk data cards 
is performed in subroutine IFS1P of the IFP module, but the detailed processing of these 
cards is done in subroutine EMGCNG of the EMG module. Besides checking for various 
errors in the CNGRNT data, the EMGCNG routine sets up a table of congruent element 
IDs in open core. This important table forms the basis for handling congruent elements 
subsequently in subroutines EMGPRO and EMGOUT of the EMG module. It is, therefore, 
useful to know the manner in which this table is set up. (The detailed manner in which 
congruent elements are handled in the EMG module can be ascertained from the source 
code and from Reference 3.) 
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The congruent element table consists of a pair of words for each element (primary 
or secondary) specified on a CNGRNT card. The first word of this pair contains the 
user-specified ID of the element. The second word of the pair indicates whether the 
element identified by the first word is the primary ID of that congruent set or is a 
secondary ID of that set. In the case of a secondary ID, the second word of the pair is 
a positive integer specifying the open core address of its primary ID. In the case of a 
primary ID, the second word of the pair is either zero or a negative integer. Initially, 
this word is set to zero in the case of all primary IDs. When the element matrices for 
the first element (which is also the lowest numbered element) in a congruent set are 
computed, the zero in the corresponding word is changed to the negated open core 
address of the element matrix (or dictionary) data. 

The second word of any pair of words in the congruent element table thus contains 
very important information. If it is a positive integer, then it is a pointer to the 
primary ID of that congruent set. If it is zero, then the corresponding ID in the first 
word is the primary ID and the element matrices have not yet been computed for that 
set. If it is a negative integer, then the corresponding ID in the first word is the 
primary ID and the element matrices have been computed for that set and can be obtained 
from the open core address information represented by that negative integer. The table 
as set up by the EMGCNG routine is sorted on the element IDs in the first words of the 
word pairs. An example of a congruent element table is shown in Table 1. 

The unique design of the congruent element table allows for a very efficient pro- 
cessing of the congruent data by the EMG module. It should be noted that the EMG 
module will never mix element matrices for different element types. The EMG module 
processes each element type one after another. When it completes the processing of an 
element type, the negated open core addresses in the congruent element table (if any) 
are replaced by zeroes. Thus, when the processing of the next element type starts, 
the congruent element table (if any) has no history or evidence of the processing of the 
previous element type. Note also that the design of the table permits the specification 
of non-existent element IDs in the CNGRNT data. 


EXAMPLES OF CONGRUENT FEATURE USAGE 


There are 82 demonstration problems in Level 17.5 of NASTRAN. The congruent 
feature is employed in fifteen (15) of these problems. A comparison of the EMG module 
CPU times (on IBM S/360-95 computer) for these problems with and without the 
congruent feature is presented in Table 2. The savings resulting from the use of the 
congruent capability are quite apparent from this table. The most dramatic savings are 
obtained in NASTRAN Demonstration Problem Nos. 3-1-2 and 8-1-2 (UMF Problem ID 
Nos. 30120 and 80120, respectively) in which the EMG module CPU times are reduced 
by more than 99 percent. 
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SUMMARY 


The congruent feature in NASTRAN is explained and the software design charac- 
teristics affecting its usage and the factors affecting its efficiency are discussed. The 
details pertaining to the software design of the capability are presented. Examples 
illustrating the usage of the feature are considered. The results of the paper clearly 
demonstrate the role of the congruent feature in increasing computational efficiencies 
and its applicability to large-size problems. 
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APPENDIX 


Input Data Card CNGRNT Identical Elements Indicator 

Description : Designates secondary element(s) identical to a primary element. 


Format and Example : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

CNGRNT 

PRID 

SECIDI 

SEC 1 02 

SECID3 

SECID4 

SEC I 05 

SECID6 

SECID7 

abc 

CNGRNT 

11 

2 

17 

34 

35 

36 





+bc 

SECID8 

SECID9 


-etc.- 

















A1 ternate Form 


CNGRNT 

PRID 

SECIDI 

"THRU" 

SECID2 






CNGRNT 

7 

10 

THRU 

55 







Field Contents 

PRID Identification number of the primary element (not necessarily the owest number) 

SECIDi Identification number(s) of secondary element(s) whose matrices will be identical 

(or congruent) to those of the primary element. 


Remarks : 1. Orientation, geometry, etc. must be truly identical such that the same stiffness, 

mass and damping matrices are generated in the global coordinate system. 

2. This feature is automatically used by the INPUT module. 

3. The CNGRNT feature cannot be used when an AXIC card is present in the bulk data deck. 

4. An element that has been listed as a primary ID on a CNGRNT card cannot be listed as 
a secondary ID on another CNGRNT card. However, if the element is listed as a 
secondary ID on the same card, then such secondary IDs are ignored. 

5. The same secondary IDs cannot be listed as congruent to two or more different 
primary IDs. 

6. Redundant specifications on CNGRNT cards are ignored. 

7. The stiffness, mass and damping matrices are actually calculated for the lowest 
numbered element in the congruent set (even though this element may not be the 
primary ID). 
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TABLE 1. - EXAMPLE OF A CONGRUENT ELEMENT TABLE 


Open Core 
Location 

Table 
Column 1 

Table 
Column 2 

Z (53) 

1 

54 

Z (54) 

3 

0 

Z (55) 

7 

59 

Z (56) 

12 

54 

Z (57) 

15 

59 

Z (58) 

18 

54 

Z (59) 

36 

0 

Z (60) 

40 

59 

Z (61) 

69 

54 


Note : The above table represents the congruent element table as initially set up 

by EMGCNG routine resulting from the processing of two CNGRNT bulk data 
cards — one with a primary ID of 3 and secondary IDs of 1, 12, 18 and 69 and 
the other with a primary ID of 36 and secondary IDs of 7, 15 and 40. 
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TABLE 2. - EXAMPLES OF CONGRUENT FEATURE USAGE IN LEVEL 17.5 NASTRAN DEMONSTRATION PROBLEMS 


Example 

No. 

Demonstration 
Problem No. 

UMF 
Problem 
ID No. 

Congruent Element Data 

EMG Module CPU Times (sec 6 )* 


Element 

Type 

Number of 
Elements 

Number of 
CNGRNT Sets 

Using the Congruent 
Feature (a) 

Without Using the 
Congruent F eature 

(b) 

by Using the Congruent Feature (%) 

i 

1-3-1 

10310 

QDMEM 

216 

1 

0.8 

8.3 

90. 4 

2 

1-3-2 

10320 

QDMEM1 

216 

1 

1.2 

13.5 

91.1 

3 

1-3-3 

10330 

QDMEM2 

216 

1 

1.5 

11.1 

86.5 

4 

1-8-1 

10810 

HEXA1 

40 

1 

0.1 

3.5 

97.1 

5 

1-9-1 

10910 

HEXA2 

40 

1 

0.3 

7.4 

95.9 

6 

1-11-1 

11110 

QUAD1 

50 

1 

0.2 

7.7 

97.4 

7 

1-13-1 

11310 

IHEX1 

40 

5 

2.8 

16.9 

83.4 

8 

1-13-2 

11320 

IHEX2 

2 

i 

2.7 

4.5 

40.0 

9 

3-1-1 

30110 

QUAD1 

200 

1 

0.4 

15.4 

97. 4 

10 

3-1-2 

30120 

QUAD1 

800 

1 

0.8 

130.5 

99.4 

11 

5-1-1 

50110 

TRIA1 

80 

4 

0.7 

11.7 

94.0 

12 

8-1-1 

80110 

QUAD1 

100 

1 

0.4 

5.8 

93.1 

13 

8-1-2 

80120 

QUAD1 

400 

1 

0.4 

49.1 

99 0 2 

14 

14-1-1 

140110 

QUAD2 

10 

5 

1.7 

2.3 

26.1 

15 

15-1-1 

150110 

BAR ) 

10 I 

5 ) 








{ 

( 

1.4 

5.0 

72.0 




QUAD2 ) 

20 ) 

5 J 





♦All of the above problems were run on the IBM S/360-95 computer® 




DEVELOPMENT OF THE LEARJET 28/29 WING USING NASTRAN ANALYSIS 


Robert R. Boroughs 
Gates Lear jet Corporation 


SUMMARY 

A great deal of the structural development work performed on the Learjet 
28/29 wing was accomplished using Nastran analysis. This included the basic 
sizing of primary structural members such as wing skins, wing skin splices, 
and spar caps, as well as the calculation of preliminary weight estimates 
utilizing the weight computation routine in Nastran. The eight spar redundancy 
of the Learjet wing made this task somewhat more complex and challenging than 
for the more determinate type wing structures. The discussion that follows 
describes some of the problems that were encountered and the solutions and 
methods that were used. 


INTRODUCTION 

The Learjet 28/29 wing was the most significantly different derivative 
wing both structurally and aerodynamically, since the introduction of the 
Learjet Model 23. This wing evolved from the earlier Model 35/36 wing and 
has been installed on a modified Model 25 fuselage. The most outwardly notice- 
able changes to the 28/29 wing from previous Learjet wings were at the wing 
tips. Here the two foot extension and tip tank on the Model 35/36 wing were 
replaced by a six foot extension and winglet on the 28/29 wing (see Figures 
1 & 2). The outward appearance of the 28/29 wing in the inboard section 
remained the same as the 35/36 wing, and internally this section still had 
eight spars as does the 35/36, but this was where most of the similarity ended. 
The wing skin and center line splice plate thicknesses have increased, the 
section properties of several of the spar caps have increased, and wing skin 
stringers have been extended or added. 

This same basic 28/29 wing configuration was later selected as the airfoil 
for the Learjet Model 54/55/56 series aircraft. Some growth capability was 
included in the 28/29 wing for the 50 series aircraft, but the complete detail 
structural definition was to be determined further into the 50 series project. 


BACKGROUND 

Approval to proceed with the development of the Learjet 28/29 aircraft 
was received in February of 1977, and work began on the wing structural analysis 
using Nastran that same month. The basic objectives established for the wing 
structural redesign were to obtain satisfactory margins of safety, minimize 
the impact on tooling, keep the weight increases as small as possible, and 
complete the certification on a very tight schedule. These goals were to be 
achieved while operating under constraints such as limited manpower availability 
and increasing lead times for parts and materials. Factors such as these later 
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influenced alternatives that were chosen during the course of this project. 

Initial analytical work performed on the 28/29 wing inboard section was 
accomplished using the 35/36 Nastran wing model described in NASA TMX-3428 
(Ref. 1). This model was later updated in the outboard section with the six 
foot extension and winglet attachment structure. Sizing of the structural 
members was to be accomplished by a combination of Nastran analysis, post 
processor programs, and detail stress analysis. Since the results of the 35/36 
Nastran wing model had correlated very closely with the experimental data from 
the 35/36 wing static test, the strategy for the 28/29 wing qualification was 
to perform limit load tests on a highly instrumented wing static test article, 
establish a correlation between the Nastran results and experimental data at 
this point, and qualify the ultimate load conditions by analysis using Nastran 
results for FAA certification. The advantage of this type of approach was to 
reduce the costs and lead time associated with a static test. 


MODEL DESCRIPTION 

Since the 28/29 wing was symmetrical about the aircraft center line, a 
half model was used. The Nastran wing model geometry was developed from the 
35/36 wing contours inboard of W.S. 181.10 and from the 28/29 wing contours 
outboard of W.S. 181.10. The wing surface was divided into a basic mesh which 
was defined by the intersection of the spar caps and rib caps. This was further 
subdivided in the spanwise direction by breaking these bays into equal incre- 
ments where possible. Structural members modeled included the spar caps, rib 
caps, and wing skin stringers with ROD elements, the spar webs and shear webs 
with SHEAR elements, and the wing skins and wing skin splices with QDMEM2 
elements (See Ref. 2). 

Modeling of the wing skins included the effect of sculpturing and contour- 
ing, and the lower skin reflected the stiffness of the access doors. The wing 
skin stringers were modeled by dividing the stringer areas in half and lumping 
each half as a separate element with the adjacent spar cap. This was done in 
order to simplify the modeling and conserve degrees of freedom. There were 
four different splice plates on the 28/29 wing skin. These were the upper and 
lower splices at W.S. 0.00 and the upper and lower splices at W.S. 181. 10 where 
the six foot wing extension was attached. The finite element representation 
reflected the contour and taper characteristics that had been machined into 
these members. 

The six foot extension geometry was basically an extension of the taper 
and contour of the inboard wing section. The inboard eight spars were continued 
into the six foot extension, and two additional spars were added in the trailing 
edge. This spar addition was incorporated to provide stiffness and an internal 
load path for forces developed by the winglet, since the winglet was mounted 
very near the trailing edge of the wing. The winglet attachment Structure was 
modeled from W.S. 244.10 to winglet station 6.00. This structure included spars 
five through ten and the winglet skin in this area. ROD elements were used to 
model the winglet spar caps, SHEAR elements were used to model the winglet spar 
webs, and QDMEM2 elements were used to model the winglet skin. Beyond winglet 
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station 6.00 the winglet structure was modeled strictly as a load fixture using 
Nastran BAR elements (See Fig. 3). 

Constraints for the model were applied in the spanwise direction at the 
W.S.0.00 spar caps, and in the vertical direction at the wing attachment 
fittings at spars 2, 5, 7 and 8, and in the fore-aft direction at spar 5. Ten 
basic load conditions were examined during the static analysis. These cases 
consisted of positive and negative gust loads as well as various landing loads. 


STRUCTURAL CONFIGURATION DEFINITION 

Preliminary design analysis of the existing 35/36 wing structure had 
revealed that larger section properties or higher allowables would be necessary 
to sustain the increased 28/29 loads. During the initial phase of wing redesign 
strong emphasis was placed on retention and utilization of existing tooling. 

This was done in an attempt to keep tooling costs and the parts count down and 
simplify the fabrication and assembly process. Consequently, several different 
configurations were analyzed where the wing skin was selected as the primary 
member for material addition with reinforcement of the spar caps in localized 
areas as the secondary means of material addition. This type of approach 
normally has been reserved to supersonic airfoil construction where a thin 
wing chord section eliminates many possible structural configurations (See Ref. 
3), but in this situation the constraints were more cost oriented. 

Each successive configuration examined had a thicker skin than the preced- 
ing configuration, but many of the spar caps still had unacceptably high stress 
levels. By this time the weight increases had become substantial and the impact 
of this parameter on flight performance had become a serious factor. Conse- 
quently, this approach was eliminated as an acceptable solution for obtaining 
the basic design goals. 

A new approach was then chosen where the emphasis was placed on increasing 
the spar cap areas in combination with the wing skin thickness as the means for 
developing satisfactory stress levels. Based on the wing skin studies that 
were conducted earlier, a wing skin configuration was selected for the 28/29 
wing. The thickness of these skins were slightly greater than the existing 
35/36 wing skins, but the total thickness was also considerably less than most 
of the other configurations previously examined. The material selected for 
both the upper surface and lower surface was 2014-T6. This was the same mate- 
rial that had been used on the upper surface of the 35/36 wing, but on the 

lower surface the 2014-T6 was used in place of 2024-T3. Selection of 
this material was influenced to a great extent by raw stock lead times 
in effect during that period, as well as the change in loads from the 
35/36 wing to the 28/29 wing. 

Using the basic wing skin selected from the previous studies, spar cap 
areas were increased in the regions where the margins of safety were deficient. 
This process initially concentrated on the wing section inboard of the landing 
gear rib where the stresses were the highest. When this region was improved 
to satisfactory levels, the process was expanded to the region outboard of the 
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landing gear rib, and from there on out to the winglet attachment structure. 

A localized buckling analysis which has been described in NASA TMX-3428 
(Ref. 1) was used to determine the non-linear effect of wing skin buckling on 
the spar cap stresses. This analysis generally required several iterations 
before a convergent solution was obtained. 


INTERNAL LOADS REDISTRIBUTION 


Redundancy in the 28/29 wing with the multiple spar construction has some 
very distinct advantages for fail safe capability, but this same asset makes 
the structure somewhat more difficult to analyze. Nastran finite element 
analysis has made this task more manageable and has permitted a better under- 
standing of this complex structure. As the first series of iterations on the 
inboard spar cap areas were approaching a convergent solution, there was ob- 
served a significant redistribution of internal loads from the previous con- 
figuration. As material was added to the critical sections, there appeared to 
be a significant redistribution of internal loads from the less critical areas 
into the more critical areas (See Fig. 4). Although stress levels decreased 
in the critical areas, these levels did not decrease linearly with the increase 
in spar cap area, and the stress levels in the non-critical regions also de- 
creased at the same time. This reduction in stress level in the non-critical 
areas may have seemed to indicate that material could have been removed from 
these areas to help reduce weight, but there was obviously another factor to 
consider. Further area reductions in the non-critical regions would have 
increased the stresses in the critical regions further, and created a need for 
more material additions in those regions. This redistribution of internal 
loads into a few key structural members also raised serious questions as to 
whether an effective and efficient fail safe qualification could be used for 
a structure defined in this manner. 


As a result of the concern for maintaining an effective and efficient fail 
safe capability for the 28/29 configuration, a different approach was selected 
for establishment of the spar cap section properties. This new approach 
emphasized maintaining an internal load and stress level balance across the 
chord section of the wing. This was accomplished by initially assigning equal 
areas to each of the spar caps from the leading edge to the trailing edge, 
and increasing each spar cap area by an equal increment for each iteration 
until a satisfactory margin of safety was achieved for the critical member. 

This worked out quite well and the weight penalties were not quite as severe 
as was seen in the first approach (See Fig. 5). After acceptable margins of 
safety were achieved in the critical spars, area reductions were then made 
on some of the less critical spars. These spars were generally located 
near the leading edge of the wing and the trailing edge of the wing. These 
spars were not generally highly loaded, and the reduction in these spar cap 
areas had little impact on the other spars. Further weight reductions were 
made by tapering the spar caps in the spanwise direction. 


A new material was. also chosen for the spar caps inboard of the landing 
gear rib in order to obtain higher allowables that were more compatible with 
the wing skin allowables and to also help reduce the weight of these members. 
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This material was7Q75-T73 extrusion, and in addition to the improved allowable 
values this material also had improved stress corrosion resistant properties 
over some of the other 7000 series aluminum alloys. 

Improvement of the upper surface capability outboard of the landing gear 
rib was another area which received considerable attention. On previous 
Learjets the spar caps in this region had been constructed from bent up sheet 
metal channels which tapered in thickness from the inboard end to the outboard 
end as opposed to the extruded caps attached to shear webs in the inboard 
section. All of the spar caps in the wing section outboard of W.S. 53 were of 
nearly equal areas; thus maintaining a fairly even internal loads distribution. 
The margins of safety for these members were generally not as deficient as the 
inboard spar caps, and other methods were used to correct the low margins than 
were used in the inboard region. 

The buckling analysis that was mentioned earlier revealed that there were 
a number of panels that indicated advanced stages of buckling. To help 
relieve the spar cap stresses existing stringers that were used on the 35/36 
wing were extended further into the outboard sections, and stringers were added 
to some bays where no stringers had existed previously. This not only caused 
the skins to carry more load in compression, but also added more basic area 
very near the outside fiber to help reduce the bending stresses on the spar 
caps. In some areas the use of wing skin stringers was not sufficient to 
obtain satisfactory stress levels and local reinforcements were added. 


PRELIMINARY WEIGHT ESTIMATES 

Preliminary weight estimates were arrived at with the aid of the weight 
calculation routine in Nastran. Two PARAM cards were inserted into the Bulk 
Data deck. The first PARAM card called out the GRDPNT option, and the second 
card called out the WTMASS feature. Accordingly, density factors were added 
to all of the material cards. A model 35/36 wing was run first to determine a 
base line weight upon which weight increases for the 28/29 wing would be 
determined. Although this routine does not include such detail factors as 
fuel sealer weights, control mechanism weights, and other miscellaneous factors, 
the preliminary weight estimates were still considered to be a reasonably 
accurate measure of weight increase over the 35/36 wing. 


DOWN BENDING ANALYTICAL QUALIFICATION 

Originally both the up bending and down bending ultimate load conditions 
were proposed to be qualified by analysis. Although there was a great deal 
of analytical work done on the up bending load condition, this load case was 
eventually qualified by static test due to the tight schedule and lack of man 
power availability that was prevalent at that time. However, the ultimate 
down bending load case was certified by analysis. The down bending load condi- 
tion was not.as critical as the up bending load condition, and the 28/29 wing 
did not require nearly as much rework for this condition as was necessary for 
the up bending condition. A comparison of the 28/29 down bending loads with 
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the 35/36 down bending loads showed that the 28/29 loads were greater than the 
35/36 loads, hut not by a large amount. Considering the material additions to 
the lower wing skin thickness, there was very good reason to expect that the 
stress levels might be very nearly the same. 

To achieve this analytical qualification a correlation was first established 
between the .35/36 wing static test strain levels and the 35/36 Nastran strain 
levels inboard of W.S. 181.10. A comparison was then made between the 35/36 
Nastran strain and the 28/29 Nastran strain levels. Im almost every location 
the 28/29 wing Nastran strain levels were less than or equal to the 35/36 
Nastran strain levels. In those areas where the 28/29 wing down bending strains 
exceeded the 35/36 down bending strains margins of safety were calculated, and 
in all cases these margins of safety were shown to be more than adequate. 

The wing structure between W.S. 181.10 and 244.10 commonly referred to as 
the six foot extension was qualified by analysis and proof load tests. 

Generally the stress levels in this section were low and had quite high 
margins of safety. The structure outboard of W.S. 244.10 consisted 
entirely of the winglet and winglet attachment structure. Due to the 
complexity of this member, certification was accomplished with a static 
test for both the up bending and down bending conditions. Typical plots 
showing the relationship between the 35/36 wing experimental data, the 
35/36 Nastran wing data, and the 28/29 Nastran wing data have been shown 
in Figures 6 thru 9. 


STATIC TEST UP BENDING CORRELATION 

There were over 400 strain gages installed on the 28/29 wing static test 
article. This was more than twice the number of gages installedon any previous 
Lear jet wing test. All the strain gages were installed on the right hand wing 
to simplify the installation and instrumentation. During the static test these 
gages were monitored on a cathode ray tube (CRT) using an interactive graphics 
program. This program provided a quick means of monitoring the status of the 
static test article and identifying areas that could become critical. Upon 
command the program would list the top 15 gages in tension and compression on 
the CRT as well as on a hard copy printer. Warnings were also issued for non- 
linear gages above a certain strain level, and individual stress versusload 
plots could be obtained within seconds for any strain gage channel. Using 
previously calculated Nastran stresses in key areas, a comparison was made 
with the appropriate strain gage channel to determine whether the test was. 
proceeding as planned. This approach proved extremely valuable in monitoring 
and controlling the static test (See Fig. 10). 

At the conclusion of the static test the strain gage results were used for 
a more detailed comparison with Nastran analytical results. The correlation 
of this experimental data with the Nastran data was generally very good for 
the majority of strain gage locations. Agreement was probably best in the 
wing section outboard of the landing gear rib where the structure was most 
uniform. Correlation in the section inboard of the. landing gear rib was good, 
but in some locations there was more noticeable deviation which appeared to 
be significantly influenced by structural cutouts and discontinuities in the 
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vicinity of the gage attachment. Figure 11 shows the strain values on the spar 
5 upper cap, and Figure 12 shows the strain values on the spar 5 lower cap. 

These values correlate quite well except in the area of the carry through 
fittings at wing station 53. These fittings were installed to maintain con- 
tinuity of the spar caps which were interupted by the landing gear rib at this 
point. 

Correlation of the strain gage values and the Nastran results were shown 
in Figure. 13 for the spar 4 upper cap and in Figure 14 for the spar 4 lower 
cap. Again agreement was generally very good except in a couple of areas. 

The first area was mentioned previously in. regard to spar 5 in the vicinity of 
the landing gear rib. The other area of some deviation was in the area adjacent 
to the.W.S. 181.10 rib. Here the wing skins were discontinuous and were spliced 
by a fingered and contoured splice plate both upper and lower. The spars were 
also discontinuous in this region, and splices were installed for spar cap 
continuity. 


RESULTS AND CONCLUSIONS 

The Model 28/29 wing development at Learjet was significantly influenced 
by Nastran analysis. Configuration development and member sizing were performed 
much more accurately and faster than could have been done previously. This 
permitted Learjet to determine the impact of changes in the wing structural 
arrangement and member section properties on stress levels, internal loads, 
and aircraft weight at a much earlier point in the wing development. 

Substantiation of the 28/29 wing was accomplished by a combination of 
testing and analysis where Nastran was the basis for much of the analytical 
work. During the static test phase of certification there were no structural 
failures in the 28/29 wing due to any of the design load conditions. Correla- 
tion of the analytical results with the experimental data was generally very 
good except in areas where there were discontinuities. The ultimate down 
bending load condition was qualified by Nastran analysis which reduced the 
cost and lead time for this segment of the certification. 
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FIGURE 3 - NASTRAN 28/29 WING MODEL 
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FIGURE 6 - STRAINS IN UPPER SPAR CAP 4 (DOWN BENDING LOADS) 
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FIGURE 10 - MODEL 28/29 WING STATIC TEST SET-UP 
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FIGURE 12 - STRAIN IN LOWER SPAR CAP 5 CUP BENDING LOADS 
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FIGURE 14 - STRAIN IN LOWER SPAR CAP 4 (UP BENDING LOADS) 



TENSIONING OF A BELT AROUND A DRUM USING MEMBRANE ELEMENT* 


C. H. S. Chen 
The BFGoodrioh Company 


SUMMARY 


A class of problem which can be solved by using the membrane element 
approach is that in which the membrane surface is either unchanged or allowed 
to undergo a given amount of transverse displacement. A problem belonging to 
the latter case is the tensioning of a belt around the drum. In this paper 
the belt tension increase due to drum edge wear or material trapped between 
the drum and the belt is investigated and some interesting results are 
obtained. In both cases the increase in belt tension Is due to the additional 
stretching of the belt resulting from the drum radius change rather than from 
the transverse defleotion of the belt. 


INTRODUCTION 


One shortcoming of the NASTRAN membrane elements is that in their for- 
mulations the coupling between bending and stretching is neglected (Ref. 1). 

In other words, the in-plane strains are Independent of the transverse dis- 
placements. As a result, the use of NASTRAN membrane elements is very 
restricted. One class of problems which permits the use of the membrane 
elements Is that in which the membrane surface is either unchanged (such as 
in a plane stress elasticity problem) or is allowed to undergo given amounts 
of transverse displacement. In either case, the applied load vectors should 
always be in the plane of the membrane element* This means that in this class 
of problems the membrane elements can not take either a normal pres stir e load 
or a concentrated load. The problem of belt tensioning around a drum can be 
classified into that particular class of problem and is used here as an 
example. 


ANALYSIS 


This paper demonstrates an application of the membrane element to the 
problem of the tensioning of a conveyer belt which wraps around a drum. The 
conveyer belt to be considered has one row of longitudinal wire cable 
reinforcement placed in a thin sheet of rubber. Conventionally the wire 
cables are equally spaced across the width of the belt. The dimensions of the 
belt oross-seotlon and the drum diameter are such that the assumption of the 


* The author wishes to acknowledge The BFGoodrioh Company for the permission 
to publish this paper. 
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belt as a membrane structure in the conventional engineering sense can be 
justified. The belt is assumed to be orthotropic linear elastic and its 
in-plane material properties are derived using composite material theory. 

In particular the well known Halpin-Tsai equations are used by knowing the 
constituent properties and the wire volume ratio (Ref. 2) , The tensioning 
of the belt considered here comes from the counterweight which is needed to 
develop sufficient friction between the belt and drum so that the belt can 
be driven when the driving pulley turns. 

We are particularly interested in investigating the stresses developed 
in the belt under the following situations* 

(1) When the drum edges wear out, 

(2) When there is material trapped between the belt and the drum, 

(3) When both the drum edge wearing and the material trapping occur. 

In all three cases there are dangers of an increase in cable tension 
which affects the belt safety. 

The trapped material considered here is in the form of a narrow strip 
wrapping around the drum for the entire region of belt-drum contact. Figure 1 
shows the modeling of a belt using isoparametric quadralateral membrane 
elements , This belt has 180 degree contact with the drum (called 180 degree 
wrap angle) , but because of symmetry only one half of the wrap angle is 
included in the model. The drum is not modeled in the problem since it is 
considered to be rigid. The belt-drum contact region extends from grid point 
1 to grid point 7 as shown in Figure 1, The fine grids are used at the edge 
regions of the belt for convenience in examining different degree of the drum 
edge wear. The width of the edge regions is arbitrarily assumed. 

The tensioning of the belt is simulated by imposing at the straight end 
of the belt (marked by grid points 11 and 18?) a specified amount of uniform 
displacement while retaning the wrapped end at the half -wrap angle (marked by 
grid points 1 and 177) unmoved. The drum edge wearing is simulated by allow- 
ing the grid points in the belt edge regions to be unrestrained in the radial 
direction. The trapped materials between the belt and the drum is simulated 
by specifying at the appropriate grid points within the belt-drum contact a 
specified amount of radial displacement. For all three cases just mentioned, 
the loaded belt can be said to have well defined surfaces and thus the 
membrane elements can be satisfactorily applied. 


NUMERICAL RESULTS 

The belt under consideration has the following dimensional and material 
property characteristics : 
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Ratio of edge region width (2d e ) to center region width (2d 0 ) = 2d e /2d c 

“ 0.35 

Ratio of belt thickness (t) to drum radius (R) * t/R = 1/22.529 

Ratio of belt thiokness (t) to belt total width (2d e + 2d ) = t/2(d a + 2d-) 

* 1/78.488 


Ratio of longitudinal Young’s modulus (E.) to transverse Young’s 

modulus (E 2 ) = E^/22 s l60,4l? 

Ratio of the height of trapped material (h) to drum radius (R) = h/R 

“1/31. 


The specific amount of displacement required at the straight end (or far 
end) of the belt to produce a given amount of belt tension is obtained by 
interpolation between several solutions. Once the element stress (0" Q ) is 
known the cable force (F) can be calculated thus, 

F = (T e EA/E 1 

where E and A are the Young's modulus and the cross sectional area of the 
reinforcing cables respectively. 


Figure 2 shows the results of the cable force when there is no drum edge 
wear, whereas Figure 3 shows that when there is complete drum edge wear. They 
are all normalized with respect to the ideal uniform cable force (F p ) when the 
drum has no edge wear and there is no trapped material, F 0 is obtained either 
by dividing the total belt tension by the number of cables or using the 
equation after the NASTRAN run. The force profiles of Figures 2 and 3 are 
those along the line of the half-wrap angle where the peaks are always found 
to be maximum. Five locations of the trapped material are examined. These 
locations are the belt centerline and four other locations at distances of 
d_/4, d_/2, 3d 0 /4 and d c away from the belt centerline. Here d Q is the half- 
width of the belt center region. 

It can be seen from Figures 2 and 3 that, for this particular case, the 
cable force (or stress) increases by about 3.3 times due to material trapping 
alone, about 1,3 times due to drum edge wear alone and about 4 times due to 
the combination of the two. These clearly indioate the danger of material 
trapping and excessive drum edge wear. Figures 2 and 3 also reveal that the 
force or stress concentration is highly local in nature and that the peak 
values appear to be independent of the location of the trapped material. 

These imply that the interactions between the cables are not very strong and 
they decay fast. This fact suggests that in the simplified analysis the cable 
interaction can be neglected. 


CONCLUDING REMARKS 


Using the problem of tensioning of a belt as an example, the present 
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paper has pointed out a class of problems whioh can be adequately solved by 
using NASTRAN membrane elements* This class of problems are such that the 
middle surface of the membrane is either undeformed or undergoes a specified 
amount of defleotion. We have studied the effects of material trapping and 
drum edge wear on the cable tension in a belt. It should be noted here that 
the material trapped between the drum and belt is of a special form that has 
the equvalent effect of increasing the drum radius. The increase of the belt 
tension really comes from the additional stretching of the belt resulting from 
the growth of the drum radius rather than due to the transverse defleotion of 
the membrane. 
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THERMAL STRESS ANALYSIS OF SYMMETRIC SHELLS 


SUBJECTED TO ASYMMETRIC THERMAL LOADS* 

Gordon R. Negaard 
ANAMET Laboratories, Incorporated 

SUMMARY 

This paper presents a study of the performance of the NASTRAN level 16.0 
axisymmetric solid elements when subjected to both symmetric and asymmetric 
thermal loading. A ceramic radome was modeled using both the CTRAPRG and 
the CTRAPAX elements. The thermal loading applied contained severe gradients 
through the thickness of the shell. Both elements were found to be more 
sensitive to the effect of the thermal gradient than to the aspect ratio of 
the elements. Analysis using the CTRAPAX element predicted much higher ther- 
mal stresses than the analysis using the CTRAPRG element, prompting studies of 
models for which theoretical solutions could be calculated. It was found that 
the CTRAPRG element solutions were satisfactory, but that the CTRAPAX element 
was very geometry dependent. This element produced erroneous results if the 
geometry was allowed to vary from a rectangular cross-section. The most sat- 
isfactory solution found for this type of problem was to model a small segment 
of a symmetric structure with isoparametric solid elements and apply the 
cyclic symmetry option in NASTRAN. 

INTRODUCTION 

Two recent studies have been conducted to determine stresses in ceramic 
radomes due to asymmetric thermal loadings. Transient thermal loads in both 
studies produced much sharper temperature gradients through the thickness of 
the shell than along the surface. For this reason, shell elements could not 
be used and it was necessary to use a formulation capable of modeling a three- 
dimensional temperature distribution. In the first study, four layers of 
NASTRAN level 16.0 twenty node isoparametric bricks (CIHEX2) of various thick- 
nesses were used to construct a radome model (Fig. 1). The thermal loading 
simulated a threat level laser irradiation. The results were found to be very 
dependent upon matching the nodal spacing to the temperature distribution and 
a problem size limitation was reached where economics prohibited creating a 
finer model which would be less sensitive to the temperature gradient. 

In the second study, the thermal loading simulated both axisymmetric and 
non-axisymmetric aerodynamic heating. The structure was modeled with a 
CTRAPAX axisymmetric element capable of handling both loading cases. The nodal 
point temperatures for an axisymmetric thermal load case are shown in Figure 2, 


*The work described herein was performed by the Aerospace Structures Informa- 
tion & Analysis Center (ASIAC) at the Air Force Flight Dynamics Laboratory, 
WPAFB, Ohio under Air Force Contract F33615-77-C-3046 . 
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which also illustrates the modeling of the radome nose-tip. The grid illus- 
trated represents the third iteration of the mesh size. The initial grid had 
the same spacing through the thickness but was several times coarser in the 
axial direction. This resulted in temperature gradients between nodal points 
of as much as 600°R. The NASTRAN results predicted unrealistically high 
stresses which were at first thought to be a function of the temperature gra- 
dient. The grid shown in Figure 2 reduced these gradients to less than 300°R, 
however, the stress levels were still not believable. A switch to CTRAPRG 
elements produced maximum compressive stresses of about 21000 psi, which 
agreed well with a SAAS III analysis and indicated that the CTRAPAX elements 

were indeed predicting erroneous stresses. These cases had been conducted 
with temperature dependent material properties so both elements were run with 
the temperature dependence removed, however this changed the results only 
slightly, eliminating this also as a possible source of the error. 

At this point, the reason for the variations in the stresses predicted by 
the two elements were unknown. A preliminary study on a hollow cylinder had 
shown almost identical answers for both a linearly and logarithmically vary- 
ing radial temperature distribution for the two elements. Various ways of 
modeling the radome with triangles and quadrilateral elements were investigated 
to determine if the problem was a function of modeling techniques. This did 
not appear to be so since all combinations of the AX elements produced similar 
stresses and the RG elements likewise produced a set of similar stresses. 
Figures 3 and 4 compare hoop stresses for the two elements and show that the 
CTRAPAX element predicts a ridiculously high stress of more then 600,000 psi 
in the same area. Since only the CTRAPAX and CTRIAAX elements could handle 
asymmetric loading, it was necessary to determine the reliability of these 
elements before continuing with the asymmetric aerodynamic heating case. 

SYMBOLS 

“ = thermal expansion coefficient 

V = poisson's ratio 
E = young ’ s modulus 
T = temperature field 
a = stress 
t = thickness 
b = radius 
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DISC ANALYSIS 


In order to determine why the CTRAPAX and CTRAPRG elements produced 
differing results for the radome model, a simplq disc model restrained at the 
outer circumference and which had known theoretical solutions was chosen (ref. 
1) . Both a quadratically varying radial temperature and a linearly varying 
axial temperature could be applied as shown in Table 1. The model and cross- 
sections of the axisymmetric model using rectangular elements are shown in 
Figure 7. For this model, both types of elements produced exact theoretical 
answers for aspect ratios varying from 1.0 to 10.0 for temperature independent 
material properties for combinations of radially and axially varying temper- 
atures. This eliminated aspect ratios and two-dimensional temperature gra- 
dients as being responsible for the differing results in the radome. There 
was no theoretical solution for temperature dependent material properties, 
however, both elements still produced identical answers. This left geometry 
of the elements as the only likely remaining candidate for the source of 
error. 

The geometry of the elements was changed so that all elements had at 
least one skewed side instead of being rectangular (Fig. 7). This model was 
then run with only the radially varying temperature distribution. These runs 
produced results that definitely proved that the CTRAPAX elements produced 
incorrect results when a non-rectangular cross-section is used. For the 
temperature independent results, the CTRAPRG elements produced results which 
matched the theoretical solution exactly, while the CTRAPAX elements gave 
radial and hoop stresses sixty to one hundred per cent too high. Even worse, 
these elements predicted axial stresses almost as high as the axial and radial 
stresses while the CTRAPRG results agreed with the theoretical solution of 
zero stress. The temperature dependent material runs predicted stresses that 
followed the same pattern but of course could not be compared to a theoretical 
solution. These results are shown in Figures 8 through 11. 

This analysis showed that the CTRAPRG element appeared to produce 
reliable results while confirming that the CTRAPAX element could not be trust- 
ed in a model requiring the use of non-rectangular element shapes, essentially 
ruling out the use of the CTRAPAX element in a model simulating a radome 
shape. Since only the CTRAPAX and CTRIAAX elements can be used for axisym- 
metric models subjected to non-axisymmetric loads, it was necessary to look 
for an alternate way of solving the asymmetric aerodynamic heating problem. 

RING ANALYSIS 

The cyclic symmetry option in NASTRAN was examined to determine if better 
results for symmetric structures subjected to an asymmetric thermal load 
could be obtained. A ring subjected to a temperature distribution of the 
form T=T 0 (R^) cos(nG) was chosen because theoretical solutions could be 
obtained (ref. 2). The model of the ring is shown in Figure 12. Both thirty 
degree and ten degree wedge shapes were examined, requiring twelve and thirty- 
six cases respectively when running cyclic symmetry. The model cross-section 
was deliberately made as similar to the previous disc analysis as possible 
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including the use of skewed elements shapes exactly as used in the previous 
section. Table 2 describes the formulation of the temperatures used as input 
and the resulting hoop and radial stresses to be expected. The axial stress 
should be identically zero. The actual stresses obtained at several radii 
are shown in Tables 3 and 4 for both ten and thirty degree wedges with skewed 
and non-skewed geometry. It can be seen that skewness had little effect ex- 
cept for the radial stress in the outermost elements in the ten degree wedge. 
The axial stresses tended to be less than ten percent of the lower of the 
radial or hoop stress except at the outer fiber. It was discovered that the 
axial stresses could be made smaller by making the ring thinner, thus approach- 
ing a state of plane stress more closely. Selected plots of hoop stress are 
shown as Figures 13, 14 and 15. It can be seen in these figures that as the 
wedge becomes narrower, it appears to approach the theoretical solution as a 
limit. 


CONCLUDING REMARKS 

In order to make a comparison of the computer costs of cyclic symmetry 
against the use of the axisymmetric elements, a ring model with the same 
geometry using the CTRAPAX elements with non-skewed geometry and temperature 
input at every fifteen degrees as shown in Table 1 was examined. This pro- 
duced almost exact theoretical answers, remembering that the CTRAPAX element 
required that the element shapes be rectangular while the cyclic symmetry 
option did not have this limitation. The following is a comparison of the 
running time on the CDC CYBER 74, using level 16.0 NASTRAN with 32 elements 
in each model. 


CYCLIC SYMMETRY TECHNIQUE AXISYMMETRIC TECHNIQUE 


(30° Wedge) (10° Wedge) 


(15° Increments) 


CM (octal) 165,000 
CP (sec) 445 
10 (sec) 333 


170,000 

1,200 

869 


250,000 

2,681 

307 


These results indicate that the cyclic symmetry option in NASTRAN is 
better suited to the solution of a general axisymmetric problem under asymmet- 
ric loading than the axisymmetric technique. A practical upper limit to the 
size problem that can be solved with cyclic .symmetry remains to be determined. 
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TABLE 1. PARAMETERS USED IN DISC MODEL 


DISC PROPERTIES 

b (OUTER RADIUS = 10.0 INCHES 
t (THICKNESS) MAX = 0.4 INCH 
t (THICKNESS) MIN = 0.04 INCH 
E = 10 7 PSI 

u = % 

a = 10' 6 IN/IN/°R 


RADIAL TEMPERATURE AND STRESS VARIATION 

T (r) = To + (Ti - To) [1 - JL ] 

b 2 

a rr (r) = -*Ea (Ti-To)[-^.-£] 
« = - » 1 Ea fTi-To)ti^.-^ ] 
(5 11 (r) = o 

AXIAL TEMPERATURE A ND STRESS VARIATION 
T (z) = 1000 z 

0 rr (z) = V 2> = 

(z) s 0 


Ea T(z) 
2 (1-u) 


TABLE 2. PARAMETERS USED IN RING MODEL 


RING PROPERTIES TEMPERATURE AND STRESS VARIATION 


b (OUTER RADIUS) = 1.0 INCH 

T (r,6») = To (■ 

-|-) 2 cos 9 
b 


a (INNER RADIUS) = 0.2 INCH 

o rr (r,0) = E 

a To B k cos 9 

t (THICKNESS) = 0.08 INCH 

= E 

a To D k cos 9 

E = 10 7 PSI 





r 

Bk 


u s 

■— 

— 



.2 

0.0 

.1202 

0£ 

0 

CD 

t 

II 

a 

.3 

.03258 

.1177 


.4 

.04408 

.1163 


.5 

.04808 

.1024 


.6 

.04696 

.07364 


.7 

.04144 

.02929 


.8 

.03170 

- .03088 


.9 

.01788 

-.1070 


1.0 

0.0 

-.1990 


TABLE 3. STRESSES IN 30° WEDGE MODEL AS A 
FUNCTION OF ANGLE AND RADIUS 


Ring with Ring with 

Unskewed Elements Skewed Elements 


RADIUS, 

(inches) 

, THETA , 
(Degrees) 

RADIAL 

HOOP 

AXIAL 

RADIAL 

HOOP' 

AXIAL 

,25 

0 

41 

187 

-4 

42 

187 

-3 


30.0 

35 

162 

-4 

36 

162 

-2 


60.0 

20 

93 

-2 

21 

94 

-2 


90.0 

0 

0 

0 

0 

0 

0 

.55 

0 

87 

158 

8 

87 

158 

9 


30.0 

76 

137 

7 

76 

137 

8 


60.0 

43 

79 

4 

44 

79 

5 


90.0 

0 

0 

0 

0 

0 

0 

.95 

0 

26 

-221 

43 

26 

-221 

43 


30.0 

23 

-192 

37 

22 

-192 

37 


60.0 

13 

-111 

22 

13 

-111 

21 


90.0 

0 

0 

0 

0 

0 

0 



TABLE A. STRESSES IN 10° WEDGE MODEL AS A 
FUNCTION OF ANGLE AND RADIUS 


Ring With Ring With 

Unskewed Elements Skewed Elements 


RADIUS 

(inches) 

THETA 

(Degrees) 

RADIAL 

HOOP 

AXIAL 

RADIAL 

HOOP 

AXIAL 

.25 

0 

46 

224 

-8 

50 

224 

-5 


30 

40 

194 

-7 

43 

194 

-5 


60 

23 

112 

-4 

25 

112 

-3 


90 

0 

0 

0 

0 

0 

0 

.55 

0 

101 

178 

0 


180 

2 


30 

87 

154 

0 

89 

156 

1 


60 

50 

89 

0 

51 

90 

1 


90 

0 

0 

0 

0 

0 

0 

.95 

0 

28 

-276 

33 

18 

-279 

33 


30 

25 

-239 

28 

16 

-242 

29 


60 

14 

-138 

16 

9 

-140 

17 


90 

0 

0 

0 

0 

0 

0 
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FIGURE 8. HOOP STRESS FOR SKEWED ELEMENTS 
(TEMPERATURE INDEPENDENT MATERIALS) 




FIGURE 9. AXIAL STRESS FOR SKEWED ELEMENTS 
(TEMPERATURE INDEPENDENT MATERIAL) 
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FIGURE 10. HOOP STRESS FOR SKEWED ELEMENTS 

(TEMPERATURE DEPENDENT MATERIAL) 
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FIGURE 11. AXIAL STRESS FOR SKEWED ELEMENTS 
(TEMPERATURE DEPENDENT MATERIALS) 
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FIGURE 12. AXISYMMETRIC RING MODEL SUBJECTED 
TO ASYMMETRIC HEATING 
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STRESS CONCENTRATIONS IN SCREW THREADS 


G. Peter O'Hara 

US Army Armament Research and Development Command 
Benet Weapons Laboratory, LCWSL 
Watervliet Arsenal, Watervliet, NY 12189 


SUMMARY 


The concept of stress concentration in screw threads is defined as a ratio 
of maximum fillet stress normalized to shear transfer rate. The data is 
presented as a plot of fillet stress vs, radial stress for a particular thread 
form. The Heywood equation is used to generate the basic plots and NASTRAN is 
used to extend the analysis to the case both where flanks of an individual thread 
tooth are in contact and the case where a finite axial stress is superimposed. 


INTRODUCTION 


The concept that stress or flow lines concentrate around various structural 
discontinuities is very old and has been the subject of many thousands of books 
and technical papers. It is convenient to express this concept in terms of a 
stress concentration factor (K) using the simple equation: 


Where K is a ratio between the maximum stress and some nominal stress, the single 
book "Stress Concentration Factors" by R. E. Peterson (ref. 1) is a compilation 
of the work included in some 378 references. The bulk of this work is contained 
on graphs which are plots of K vs. some geometry factor and most use a family of 
curves to show the effect of some other geometry factor. These plots provide 
both useful numeric information and a quick visual picture of the structural 
response. 

The concept of stress concentration in screw threads is rather elusive and 
in fact there is little work done on stresses in threads. R. B. Heywood (ref. 2) 
published an empirical equation for the maximum fillet stress which was used in 
the work of Weigle, Lasselle and Purtell (ref. 3) as a guide in trying to improve 
fatigue life of cannon breech rings. Later this author demonstrated that the 
Heywood equation would give accurate numeric data when the boundary conditions 
were closely controlled (ref. 4). 

However, most work with screw threads seems to be done for specific cases 
such as the fine work of M. Heyinyi (ref. 5) who investigated bolt shank and 
nut design in Witworth threaded bolts. This type of analysis using three 
dimensional photoelasticity was also used by W. F. Franz (ref. 6) and J. D. 
Chalupnik (ref. 7) . A further attempt at optimizing a thread form was done by 
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R. L. Marino and W. F. Riley (ref. 8). 


In all of these works the calculated stress concentration factor is differ- 
ent for each thread in the system. It would seem that if the stress concentra- 
tion factor is properly defined something should be a constant for all threads 
of a specific shape. In his original paper, Heywood demonstrated part of the 
problem. The stress in the fillet is the result of two factors. First, is the 
stress due to the load on the individual thread tooth and second, is the stress 
due to the general stress field or the axial stress (cr a ) near the thread fillet. 
In this paper I will add effect due to friction and normalize all stresses to 
the average shear transfer rate (x R ) . 

When the friction force and the force due to the "wedge" effect of the 
loaded flank of the thread are combined a radial (normal) force is produced 
which can be averaged into the radial stress (cr r ) . The fillet stress (Op) can 
be expressed as the sum of two functions. 

c F / T r = c F = G 1 (a,g,R,e,a r ) + G 2 (a,£,R,0,cr a ) 

In the above equation the first function (G^) is the relation between 
fillet stress and the load on the individual thread tooth. The second function 
(G 2 ) is the factor due to the general stress field. Alpha (a), beta (3) and R 
are the thread geometry factors. The angle (0) is in both functions because 
they do not maximize at the same position in the fillet. In this paper the 
shear transfer rate is defined as the net load supported by the thread divided 
by the area at the pitch line. The direction of the net load is parallel to 
the pitch line and in the analysis this component of the force will be unity. 

The radial stress (a ) and axial stress (a ) are normalized to the shear transfer 
rate. 


The above discussion relates to a normal screw thread problem where only 
one flank of a particular thread contacts one flank of a mating thread. In 
some structures the relative displacement in the radial direction across the 
threaded connection is such that the radial gap in the threads is closed and 
both flanks of each thread carry load. Under these conditions the radial 
component of the loads add together to produce high negative or compressive 
radial stress across the joint. The axial loads oppose each other and the 
pressure on the primary flank must become very high to overcome the secondary 
flank load. This is not a common condition; however, it may become very impor- 
tant in the cannon breech- to- tube connection. 


THREAD GEOMETRY 


The thread geometry parameters are shown in Figure 1 and in this report 
all linear dimensions will be normalized to pitch (P) . The primary geometry 
parameters are the primary flank angle (a), the secondary flank (3) and the root 
radius R. These, in conjunction with the pitch space (PI) , define the basic 
thread geometry. Other factors are required to insure a practical thread which 
will fit together. The addendum (AD) and dedundum (DD) dimensions sum to the 
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total height (HT) . The tip radius (RT) eliminates a sharp corner and helps to 
support the bearing surface. The root flat (FLAT) is often used to make up for 
the root radius tolerance. The bearing height (Z) is used to calculate the 
average bearing stress and the shear length (S) is used to calculate the maximum 
shear-out failure load. 

This complicated system is simplified by the fact that we must deal with 
a small set of standard thread forms. In this report detailed analysis has been 
done on the British Standard Buttress thread and Heywood analysis has also been 
done on the controlled root bolt thread or "J" thread and the Watervliet Arsenal 
Buttress used on cannon breeches. The nominal dimensions for these threads are 
shown in Table I. 


TABLE I. THREAD GEOMETRY DEFINITION 



British 

Buttress 

Watervliet 

Buttress 

30 "V" 
(Rolled) 

a 

7° 

20° 

30° 

6 

45° 

45° 

30° 

R 

0.1205 

0.1333 

0.1804 

Pi 

0.500 

0.5276 

0.500 

HT 

0.5059 

0.4787 

0.6077 

RT 

0.00 

0.0480 

0.1083 

FLAT 

0.0 

0.0 

0.0 


LOADING PARAMETERS 


The Heywood load parameters are also shown in Figure 1. They are a point 
force (W) applied at some position (b) in the loaded flank with a friction 
angle (y) . This scheme can be repeated many times on the loaded flank to 
produce some load distribution curve. The following loading assumptions are 
made: 

1. The total load vector parallel to the datum line is unity. 

2. The load distribution is uniform. 

3. Friction does not vary along the flank. 

The first assumption given allows the normalization of stresses to shear 
transfer rate and the other two establish a simple loading case. 

Under conditions of high radial compressive load, it is possible for 
threads to be pushed together until both flanks contact. This condition will 
be discussed later. Under normal conditions only the primary flanks contact on 
the thread and the radial stress become a function of the flank angle a and the 
friction angle y: 

a Y = tan (a-y) 
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Note that friction becomes a signed variable depending on the relative displace- 
ment of the two components of the structure. 

In the above discussion the general field or axial stress is assumed to be 
zero. In the NASTRAN finite element analysis the axial stress was simulated by 
the use of a constraint subcase in which the relative axial displacement 
between the two radial boundaries was fixed by the use of scalar points and 
multipoint constraint equations. The radial displacements on these planes were 
made equal for congruent points. The radial displacement of the inner axial 
boundary was set equal to the Poisson contraction of a solid bar. 


HEYWOOD ANALYSIS 


The Heywood equation is shown in Figure 2. This is a semi-emperical equa- 
tion that was fit to a large body of photoelastic data. It calculates maximum 
fillet stress for a point load on the primary flank of a thread using a specific 
friction angle. In order to simulate a uniform load distribution, the results 
are averaged for seven different "b" values evenly distributed over the flank. 
The process has been programmed info a program called HEY40. The calculations 
have been done for many standard thread forms, and the three reported in Figure 
3, have been defined in Table I. This plot of fillet stresses plotted against 
radial stress will be referred to as the "thread characteristic curve". This 
curve covers a friction angle range of -45° to 45° or a coefficient of friction 
range of -1.0 to 1.0. 

In Heywood' s photoelastic experiments he was careful to transfer the load 
supported by the threads profiles in a shear mode to make the axial stress as 
small as possible. This process limited his equation to the case where axial 
stress is equal to zero. 


NASTRAN FINITE ELEMENT ANALYSIS 


The finite element work was done for three reasons: (1) to verify the 

Heywood analysis; (2) to examine the two-flank problem; and (3) to include a 
finite axial stress. The grid for the British Standard Buttress is shown in 
Figure 4. It contains 216 triangular ring elements (CTRIARG) and 133 grid 
points. The run required five basic loading subcases plus fourteen subcase 
combinations for each value of axial stress. These fourteen subcases cover both 
1-flank and 2-flank contact over a range coefficient of friction of -1.0 to 1.0 
in increments of 0.25. 

The grid was generated using IGFES(ref. 9) and following that, force sets 
were calculated to apply uniform pressure and uniform shear loads on both flanks 
of the thread and a displacement was calculated for a nominal 1.0 psi axial load 
on the grid. Two different constraint conditions were required to complete the 
boundary conditions for a single thread taken from a long series of threads. 

For loads on the thread the inner boundary points were fixed in both radial and 
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axial directions and similar points on the two radial boundaries were con- 
strained to equal displacements. In this way the net load was taken out as 
shear load on the inner axial boundary and the multipoint constraint equations 
replaced adjoining material. For the axial load condition the inner axial 
boundary was constrained to the Poisson displacement in the radial direction 
and left free in the axial direction. The two radial boundaries were given 
fixed relative axial displacements and the radial displacement was made equal 
for similar grid points. This condition was set to simulate a far removed 
axial loading. 

Because the basic loads were all for a 1.0 psi uniform applied pressure or 
shear and the results were desired for a 1.0 psi shear transfer rate (calculated 
at the datum line), it became necessary to calculate the correct Subcase Sequence 
Coefficients for fourteen subcases for each of four axial stress values (or 56 
sets). Therefore a small program was generated to supply all necessary SUBCOM, 
SUBSEQ and LABEL cards for that portion of the case control deck. 

Uniform increments of 0.25 in coefficient were used from -1.0 to 1.0. If 
the friction value was positive, a single subcase combination was generated 
which would superimpose the two primary flank loads with the proper value of 
axial stress. If the friction was zero or less, a similar subcase combination 
was generated along with one where both flanks were loaded and the second radial 
stress was 1.0 greater than the initial. 

The axial load subcase produced the conventional fillet stress concen- 
tration factors of K = 2.89. This maximum stress was in an element at the 
bottom of the fillet where 0 is approaching 0°. In the cases where the load 
is applied to the thread the exact position of the stress maximum is about 45° 
up from the bottom of the fillet. Data is reported here for two cases of axial 
stress. Zero axial stress is shown in Figure 5 and axial stress of 2.0 is shown 
in Figure 6. These plots are a set of six lines with the single contact curve 
at the right. Starting at that curve is a family of five lines going to the 
left which represent two flank contacts at different values of friction from 0 
to -1.0. 


DISCUSSION 


The first thing to note is that there is excellent agreement between Heywood 
and NASTRAN over most of the range of the plots for the one case in question. 
Because of this, the Heywood relation can be used to evaluate different thread 
forms. The finite element method has allowed the expansion of the basic plot 
to the 2-flank contact problem and the addition of the axial stress. 

There are several important points that are demonstrated by this work. 

Note that a small change in friction can produce a large fillet stress variation 
in all three threads reported in the Heywood analysis. Negative friction angles 
can produce marked reductions in thread fillet stress. This effect was noted 
several years ago in an unpublished three dimensional photoelastic study where 
the model was overloaded and the threads were forced to a high negative radial 


69 



stress. In this case the fillet stresses were very low and the experiment was 
repeated. This author suspects that friction variation may be responsible for 
much of the scatter in bolt-fatigue data. 

This work was initiated becuase of the necessity of analyzing a structure 
with a long threaded connection using many small threads. In this case the 
modeling of each thread would require an excessively large data deck. There- 
fore, the threads were handled as a conventional contact problem where friction 
could take on any value and limits were applied to the radial stress. In the 
solution the contact surface was placed on the datum line of the threads and 
one or two teeth were replaced by one element space. Shear transfer rates could 
then be estimated from the shear stress data near the contact surface along with 
radial and axial stress. The fillet stresses were estimated for use in fracture 
mechanics analysis. 


CONCLUSION 


A stress concentration approach to the thread fillet stress problem has 
been defined using the shear transfer rate as the fundamental quantity. This 
stress concentration is plotted for a fixed geometry in a stress vs. stress 
plot where the stress concentration is a function of the applied radial stress. 
This process can be repeated for several values of the applied axial load. The 
effects of axial stress and applied thread loads seem to be of equal importance 
and accurate results require the analysis of both factors. 
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CONDENSING LOADED POINTS FOR 'TRANSIENTS BY S UB STRUCTURING ' 


THOMAS G. BUTLER 
BUTLER ANALYSES 

SUMMARY 


This paper is an extension to my article presented at the 7th 
NASTRAN Colloquium entitled "TRANSIENTS BY SUBSTRUCTURING". 

In order to capitalize on the economies of (a) condensing of 
omitted points during substructuring, and (b) of applying BANDIT 
at the only place where it can enter the analysis, a strategy 
needed to be worked out so that dynamically loaded points could 
also be condensed in Phase 1. This would obviate having to re- 
tain points during dynamic analysis for no other reason than to 
make them available for assigning load, and having to suffer the 
consequences of expensively large matrices. This paper describes 
a technique, vintage 17.5, for accomplishing these aims. A num- 
ber of problems arose such as transferring of condensed loading 
information from R. F. 1 to R. F. 9, and giving this loading a 
correct time history. The method has been applied to substructure 
transient solutions according to the approach as outlined at the 
7th NASTRAN Colloquium. Only DMAP statements and some manual 
transfer of data were used to implement this strategy. 


PURPOSE 


The cost of solving large order transient problems can be 
high. If the matrix order is reduced by condensation without in- 
flating the band, the cost of solving transients can be decreased. 
The cost of condensing in itself can be expensive for large ma- 
trices, so it behooves one to condense at the component level 
whenver the matrices are small and the cost is exponentially 
less. If this condensing of dynamically loaded degrees of free- 
dom could be worked out, then the rder and band of the matrices 
in the transient solution could also come under the control of 
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the analyst to conform to his needs without operational con- 
straints. The avenue is open in substructuring to do the con- 
densing of each low order component matrix for small cost with a 
net overall saving. It appears that the ideal is at hand to cut 
the cost of transients without incurring high processing costs, 

but the way is not straightforward. The plan should not be coded 

without considering the needs of banding. One is compelled to 
condense at the compoonent level when possible, because BANDIT 
can be invoked in substructuring only at the component level. 

The purpose of this enterprise is to find a way of condensing dy- 
namic loads in stubstructuring economically, so that the load 

matrices can be correctly formed when R. F. 9 is entered. 


PROBLEMS 


How can the dynamic loads be represented in Phase 1? How 
can the spatial, time amplification, and time delays be coordinated, 
if the points scheduled for condensation have loading times that 
differ from those of retained points? How can loading on the re- 
tained points be merged with that resulting from condensing opera- 
tions on the omitted points? Can all these objectives be met with- 
out burdensome labor or cost? The answers to these questions will 
be taken up in the section entitled IMPLEMENTATION. The reason 
for having to pose the first question is that Phase 1 of fully 
automated substructuring is operational for statics and eigen- 
values only and not dynamics. Consequently, some property of the 
dynamic load has to be isolated which behaves in a way that can be 
treated in statics. The answer to the second question seems to 
be even more elusive than the first. It appears like an attempt 
to attach a time variation to something which has disappeared from 
the problem. The merging difficulty posed in the third question 
seems to be one of retaining some sort of identity of the con- 
densed loads so that the relationships to the loads on uncon- 
densed points can be maintained. The bulk of this report is spent 
on how these questins were answered. 


IMPLEMENTATION 


Consider the composition of the dynamic load in its simplest 

form 


P(t,x,y,z) = A(x,y , z) x T(t-p) 



The component, A(x,y,z), is the spatial specification of the load- 
ing. The data for this part is entered via DAREA cards. In that 
the information contained on DAREA involves position and magnitude, 
it is similar to stati loads. If one were to set up the TL0AD1 
cards for the dynamic loads without regard for any intention to 
condense loaded points, it would be easy to identify those spatial 
sets (i.e. DAREA sets) that are loaded synchronously; i.e. those 
having a common time amplification defined by a single TABLEDl set 
and a single DELAY set. At completion of defining the total dynamic 
load there will be as many synchronous sets as there are TL0AD1 
sets. There are, at most, the same number of DAREA sets as there are 
TLOADl sets — in general there are fewer DAREA sets. It is now a 
matter of comparing the omitted degrees of freedom with the DAREA 
sets to find their intersection. All those DAREA sets having omit- 
ted degrees of freedom are isolated. For purposes of illustration, 
concentrate on just one substructure and assume there are three 
such DAREA sets in that one substructure and designate them as 
Q,R & S. Use lower case letters q,r & s for the sets of omitted 
degrees of freedom in their respective parents. Delete q,r & s 
from Q, R & S and label these DAREA sets with their condensed d.o.f.'s 
removed: Q 1 ,R' & S'. The case in which all degrees of freedom 
in one of the DAREA sets are omitted is admissable; for instance 
q could be identical to Q and Q' could be null. With this notation 
problem behind us it will ease our discussion of the strategy. 

These individually synchronous sets (q,r & s) that are both loaded 
and omitted, will be operated on during substructuring. In effect 
q,r & s are going to be condensed as static loads in Phase I then 
delivered back into the dynamic loads as pr e-condensed DAREA sets. 

The logic follows from the fact that spatial condensation in statics 
operates exactly the same way that spatial condensation works in 
dynamics. Time varying amplification won't operate any differently 
on a condensed set that comes pr e-condensed from statics and is not 
further condensed in dynamics than it would on one that comes 
earmarked to become condensed during transient execution. Con- 
sequently, when the loads q,r & s are condensed in substructuring, 
the resulting redistribution of loads to retained points will be 
labeled q',r' & s'. The next step is to organize data from q', r' & s' 
in DAREA format and give them the same set I.D. 's respectively 
as the parent Q' ,R' & S' sets. Belonging to the same DAREA set, 
constitutes a merge of sets in Case Control management, so the 
q,r, & s that were deleted from Q,R & S will have been replaced 
with their pr e-condensed counterparts q',r' & s'. No changes 
are made to the dynamic loading data involving TLOADl, DLOAD, 
and TABLEDl entries. New DELAY cards will be organized. The 
sets of synchronously omitted degrees of freedom q,r & s will 
have to be replaced by the corresponding sets of retained degrees 
of freedom q',r' & s' without any change to the delay times. 

When more than one substructure has condensed loading, the same 
procedure as outlined above should be followed for each component. 
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So much for talk. It remains to be seen whether this idea has 
substance or is just a brave show. The DAREA entries of loading 
coefficients that would have been assigned to the omitted d.o.f.'s 
of a substructure in a solo transient run are assigned instead to 
those same omitted d.o.f.'s, as static force loads in a Phase I 
run of that substructure. Each synchronous set is organized as a 
static load set and each is scheduled in a succession of subcases. 

In terms of the symbols of the previous paragraph, the DAREA 
entries on the set q wil be arranged in a static load set for the 
first subcase, those on r for the second subcase, and those on s 
for the third subcase. Of the three solution routes available in 
automated substructuring, the choice converges on R. F. 2 (Inertia 
Relief) , because it offers options for stiffness, mass, and load 
matrix generation, while R. F. 1 omits the mass and R. F. 3 omits 
the static load. If one scans the DMAP sequences for R. F. 2 and 
compares it with the automated ALTER statements for SUBSTRUCTURE 
one sees that the load vector is condensed down to retained degrees 
of freedom in the module SSG2. It is output as the data block PL. 
NASTRAN needs to be told that all 3 matrices are wanted. In the 
SUB SCC this is indicated by selecting OPTIONS = K,M,P. Nothing 
more is essential during Phase I except to emphasize that BANDIT 
should be enlisted prior to Phase I to resequence all points ex- 
cept mating points at the interfaces between substructures. The user 
should manually assign the sequence numbers to interface mating points 
Even though all the Phase I requirements are met, we are not going 
to leave it with such an uncaring attitude. Something extra is 
going to be worked up to give a better "feel" as to how things went 
during Phase I. We are going to "gin up" an ALTER packet to see 
what the condensed DAREA components look like before being content 
with the substructure component runs. The ALTER packet will use 
the module SOFI to pick off PL, as it is delivered to the SOF, and 
then will inflate PL with zeroes in a succession of merges until it 
is G-sized. Then this matrix of "load vectors" is delivered to SDR2 
for processing the OLOAD requests in Case Control and reformating 
according to external grid sequencing. Now the OFP handles the 
printing chores. Particulars of the DMAP ALTER packet are given 
in Table I. The above outlined procedure is carried out for every 
substructure which contains omitted loaded points. 

If further condensing is decided on in Phase II, using the 
REDUCE command, the bookkeeping of the several substructures is kept 
straight inside the SOF data base. When the Phase II operations 
have prepared the final pseudo structure, the SOF will contain all 
of the information that is needed to proceed into transients, pro- 
vided that OPTIONS = K,M,P has been used for every Phase II run. 

The stiffness and mass matrices are in final form and the PVEC 
item of the SOF contains a vector for each synchronous subcase of 
coefficients for all component substructures in fully condensed 
form. Now an ALTER uses the module SOFI to take the stiffness, 
mass, and load data from the SOF and deliver them externally. Use 
0UTPUTT1 for K & M for subsequent insertion into R. F. 9. Use 
0UTPUT2 for P if a FORTRAN preprocessor is to be used. If DAREA 
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input is to be prepared manually then the PVERC matrix need only 
be printed. MATPRN will be sufficient since the DAREA data will 
be written by internal scalar index number instead of external grid 
point numbers. Details of the Phase II ALTER are presented in 
Table II. The highlights of the strategies for the Phase I and 
Phase II solutions are depicted in Tables V & VI. 


EXAMPLES AND CHECK CALCULATIONS 


A small two component substructure problem was used to 
illustrate the method of substructure condensation of dynamic load 
prior to transient execution. Component BOX is a parallelepiped 
of BAR elements with BAR diagonals on the rear end to make a con- 
nection in the middle. Component END is a two bar appendage. 

The combined pseudo structure is called ALL. Figure 1 shows P/S 
ALL with grid points numbered and loading indicated. Grid points 
2 and 4 are loaded simultaneously with load F(t) and GP ' s 6 & 8 
are loaded simultaneously with load L(t). GP's 1 & 3 are loaded 
simultaneously with load F(t) delayed an amount T. GP's 5 & 7 
are loaded simultaneously with load L(t) delayed by an amount T. 
Figure 2 has plots of loadings L(t) and F ( t) . The direction for 
both dynamic loads F(t) and L(t) is in component direction 1 as 
shown by the double shafted arrow in firgure 1. If DAREA and DLOAD 
coefficients are held at unity then the values from figure 2 will 
become the TABLED1 entries. This pilot problem was deliberately 
kept small to test only the new features and not redo all features 
of transients by substructuring. D.o.f. no. 1 was omitted from point 
6 of load L ( t) and from point 3 of load F(t) and is indicated with 
tiny circles in figure 1. 

In the spirit of keeping things simple, only displacement re- 
sponses in the uncondensed component END will be examined as in- 
dicated by open brackets in figure 1. Since displacements can be 
recovered as part of the R.F. 9 output, there will be no need for 
scheduling Phases 3 & 4. In order to test the method, a separate 
transient analysis entirely within R. F. 9 was run. The problem 
flow is diagrammed starting with a Substructure Tree in figure 3 
and flow charts of transients first with a substructure preface in 
figure 4 and then in a stand-alone mode in figure 5. 

A measure of the chore of converting the PVEC from substruc- 
turing into DAREA data for transients can be obtained by looking at 
the effects of condensing. Figure 6 shows the DAREA values for the 
synchronous set of grid points 6 and 8 as it was prepared for the 
solo transient solution. Figure 7 is a sketch of the points to 
which the omitted loads of GP 6 and GP 3 were vectored by the SMP 
module in Phase I. Figure 8 shows the results of manually prepar- 
iing the DAREA data of the condensed loading for the synchronous set 
of grid points 6 and 8 which are now organized according to in- 
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terior scalar points. It is easy to imagine how burdensome this 
manual preparation can become for a large structure. A processor 
program is an obvious necessity. Only one digit of data was used 
to illustrate the relative magnitudes of the redistribution. The 
major contributions are circumscribed. Translations are enclosed 
in boxes and rotations are enclosed in circles. The translation 
load of GP 6 is essentially reduced to a translation at GP 5 
(Scalar dof 7) and rotations about GP 8 (Scalar dof 18) and 
GP 2 (Scalar dof 35) . The box about scalar dof 13 identifies the 
uncondensed loading of GP 8. 


RESULTS 


In order to test this method of condensing dynamic loads bv 
substructuring, a separate analysis was performed in the conven- 
tional way by submitting the entire structure with loads defined 
externally and all other features retained, using R.F. 9 in a 
solo run. The quality of the substructure analysis was based on 
making comparisons of only displacements at the tip of the appen- 
dage where the amplifications would tend to be the greatest. 
Comparisons of stresses and forces in elements were omitted for 
reasons of brevity. Displacement comparisons could be made dir- 
ectly from the two transient outputs, but stresses and forces 
would have required two extra steps plus additional work in 
Phase II. Results of the two displacement time histories are 
assembled side by side for each translational and rotational com- 
ponent in Tables III & IV. Where the displacements are large, 
the outputs match in 6 digits and start to vary in the seventh 
place. That is, the axial and vertical translations and the ro- 
tations about the transverse axis compare almost exactly. But the 
small displacements don't even match up in the first digit. In 
theory all of these extremely small displacements should have been 
zero, because the load was symmetrical; however the condensation 
was unsymmetr ical. This was a very small model that would be 
highly sensitive to irregularities, whereas in a large structure 
much smoothing would take place. This unsymmetr ical condensation 
caused numerical noise to creep into both the solo and the sub- 
structure solutions, giving non-zero values 8 orders of magnitude 
less than the signifigant displacements. Noise is highly depen- 
dent on the sequence of operations and the sequences in the two 
cases under scrutiny were different, so the noise from the two 
can be expected to be different. The net opinion is that the 
results compare favorably. 

On the basis of good correspondence in the responses from 
the two approaches, we can say that the technique of spatial con- 
densation of DAREA coefficients by simulating them as static for- 
ces is satisfactory. The method of condensing dynamic loads by 
substructuring has been established. 
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RECOMMENDATIONS 


When transients are automated for substructuring, it appears 
that the elements contained in this non-automated approach could 
be used as a basis for the algorithm. The Phase I bulk data could 
include a DAREA card, called say DAREAS, which would respond to the 
user's input of condensed loading sets. A caption in the LODS 
item could append a D to every load subcase that originated from 
DAREAS input. At transient execution time the user would indi- 
cate which sets of DAREAS data would carry the same set I.D. as 
uncondensed transient DAREA data, and which sets of DAREAS data 
would be distinct. Similarly, the Phase I bulk data could in- 
clude a DELAY card, called say DELAYS, which would respond to the 
user's input of condensed loading. Then the usual preparation of 
the elements of dynamic loads would proceed as is now done in 
transients. If further condensation is desired in transients 
when preceded by substructuring, then caution would have to be 
exercised in data recovery. Instead of passing the UDVT matrix 
to Phase IV, the data recovery would have to reconstruct the 
response histories to the full order of the original Scalar 
Point compliment, before passing the displacements to Phase IV. 
These displacements are contained in data block UPV. Modifica- 
tions to the DMAP ALTER' s for R.F.9 as outlined at the 7th Col- 
loquium would be needed to allow for an option to output either 
UDVT or UPV based on the condition of whether OMIT's are present. 
The subsequent partitioning would have optional input of either 
UDVT or UPV. 

In effect, the user could be in full command of the conden- 
sation of his problem in substructuring by using this proposed 
technique. He could condense in Phase I where it is most eco- 
nomical and where he could take maximum advantage of BANDIT. 

He could condense again in Phase II where interface points 
could be removed. Lastly, he could condense in transients where 
the final trade-offs are reviewed between the transient solution 
costs and additional condensation and data recovery costs. 
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S/S - R.F. 9 


SCALARS 

1 THRU 54 

DAREA 

2,4 LOAD 
5,7 LOAD 

COND 6-SET & 8 LOAD 
COND 3-SET & 1 LOAD 

DLOAD, 

TL0AD1, TABLED1 (UNMOD) 

DELAY 

5,7 - SET (UNMOD) 

COND 3-SET & 1 TAU 

ALT 

INPUTT1 KFiTX, IfITX 
ADD K4GG 

OUTPUT 

DISP - 11 
OLOAD - 9 


Figure 4 
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SOLO 

COMPLETE TRANSIENT 


BULK OF ALL (BOX & END) 
LOAD EXTERNAL BY GP & COMP 
OMIT 

OUTPUT DISP - 11 
OLOAD - 9 


Figure 5 



DAREA 

BEFORE OMIT 
EXTERNAL SEQ. 


G.P. 6 

SID GP COMP 
68 6 1 


G.P. 8 


68 8 1 


Figure 


VALUE 

1.0 


1.0 
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Figure 7 
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DAREA 

AFTER OMIT 
INTERNAL SEQ. 


SID 

GP 

COMP 

VALUE 

GP 

COMP 

VALUE 

68 

7 

0 

0 

8 

0 

.02 

68 

9 

0 

-.02 

10 

0 

.0 

68 

11 

0 

-.1 


12 

0 

-.1 

68 

13 

0 1 

1,1 

l 

14 

0 

-.02 

68 

15 

0 

.002 

16 

0 

-.01 

68 

17 

0 

.06 

18 

0 

© 

68 

31 

0 

.1 


32 

0 

-.002 

68 

33 

0 

.02 

34 

0 

.01 

68 

35 

0 

© 

36 

0 

.06 


□translation O Rotation 


Figure 8 



PHASE I ALTER 
PACKET 


ALTER 

SOFI 

UMERGE 

UMERGE 

UMERGE 

UMERGE 

CHKPNT 

SDR2 


SAVE 

OFP 


136 




/Pn,,,,/C,N,1/C,N,BOX/C,N,PVEC $ 

PL -jftrznrt' SQP 

USET, Pn, /PLTOA/C, N, A/C, N, L/C, N, R $ 
USELPLTOA, /PAT0F/C,N,F/C,N,A/C,N,0 $ 
USET,PATOF, /PFTON/C, N, N/C, N, F/C, N, S $ 
USET,PFTON, /PGG /C,N,G/C,N,N/C,N,M $ 

PL 

PGG $ 

CASECC, CSTM, MPT, , EQEX I N, S I L, , , BGPDT, , , , , , PGG/ 
0PG1, , , , , /C, N, STAT I CS/C, N, N0S0RT2 $ 





N0S0RT2 


0PG1, , , , , //V, N, CARDNO $ 
SAVE CARDNO $ 


ENDALTER 


Table I 
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PHASE II ALTER 
PACKET 


ALTER 136 END & 

SOFI /Kn^In, Pn, , /C, N, 1/C, N, name/C, N, KMTX/C, N, MMTX/C, N, PVEC $ 

2L s&rzul < sir6c£ry' PVf6 

kl^TUMJ X P/U&fA COM4' Mi sy\ S/S 

0UTPUT1 Kn,Mn,^//C,N,-1/C,N,INP#/C,N,(C0NDKM) $ 

ft yp 
AL^rbi^Asipy^ c^~rC 

^/T f, 

MATPRN Pn„„// $ , 

c&itPnf Jrr "tH‘ 

0UTPUT2 Pn, , , , / A?, N, -1/C, N, 11/C, N, ( LODOM I T) $ 

0UTPUT2, ,,,,//C,N,-9/C,N,ll $ 

)?U£*A 0Tp DAPFA Czryyi/urYiAvUa 

A4j s&t-rti AOF7~AAA/ s&rziM;sL.f ^ * 

TAP*! firr JrreA. J PAPFA 

ENDALTER ^ 


Table II 
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COMPARISONS OF DISPLACEMENT 
RESPONSES AT END TIP (G.P.ll) 
BETWEEN COMPLETE R.F.9 & S/S R.F.9 
AT EVERY SECOND FOR 7 SECONDS 

T1 T2 


TRC 

S/S 

TRC 

S/S 

-4,000940+2 

-4.000941+2 

-3.4-5 

-7,8-5 

-7,335053+2 

-7.335056+2 

-2.2-4 

-Li-4 

6,445970+2 

6.445966+2 

-5.4-4 

-LI-4 

5.479066+3 

5.47906Z+3 

-7.3-4 

+2,6-4 

1.451452+4 

1.451452+4 

-9.9-6 

+1,0-3 

2 , 738419+4 

2 . 738421+4 

+3,2-4 

+2,1-3 

4.287670+4 

4 . 287672+4 

+9.2-5 

+1,5-3 


TRC 

T3 

S/S 

TRC 

R1 

S/S 

3,159634+2 

3.159632+2 

-1.4-6 

+1,5-5 

5,792659+2 

5,792651+2 

-1.5-5 

+6.5-5 

3.949432+1 

3 , 949254+1 

-5.7-5 

+1,6-4 

-1.145368+3 

-1.145371+3 

-1.0-4 

+3,2-4 

-2,290734+3 

-2.290735+3 

-1,0-4 

+5,4-4 

-2,317057+3 

-2.317051+3 

+2.2-6 

+7,6-4 

-1.342825+3 

-1,342805+3 

+2.2-4 

+9,6-4 


Table III 
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COMPARISONS OF DISPLACEMENT 
RESPONSES AT END TIP (G.P.ll) 


BETWEEN COMPLETE R.F.9 & S/S R.F.9 
AT EVERY SECOND FOR 7 SECONDS 


R2 R3 


TRC 

S/S 

TRC 

S/S 

-1,182928+1 

-1.18292Z+1 

-9.0-7 

-1.0-6 

-2.167783+1 

2 . 167781+1 

-5.9-6 

+2.9-6 

-1.977996+0 

-1.97799Z+0 

-1.9-5 

+1.7-5 

9.286305+1 

9.286311+1 

-1.9-5 

+9.5-5 

8.572603+1 

8.572609+1 

+9.6-6 

+8.0-5 

8.671105+1 

8.6710Z9+1 

+6.3-6 

+1.0-9 

5.025239+1 

5.025163+1 

-1.7-5 

+3.0-5 


Table IV 



PHASE I STRATEGY 
FOR DAREA CONDENSATION 

EXECC SOL 2,0 

ALTER 

SUBSCC OPTIONS = K,M,P 

CASECC OLOAD = ASET D.O.F.'s 

S/C 1 

LABEL = DAREA GP A COMP 
LOAD = 1 

S/C 2 

LABEL = DAREA GP B COMP 
LOAD = 2 


S/C N 

LABEL = DAREA GP M COMP 

BULK F0RCE1 . , . . GP A . . COMP 

F0RCE1 . , , , GP B , . COMP 
F0RCE1 . . . . GP M . . COMP 


Table V 


VALUE Rj_r 
VALUE Sj_s 


VALUE LlT 

. . VALUE Rj_R 
, . VALUE S.S 
, . VALUE tj: 
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PHASE II STRATEGY 
FOR DAREA CONDENSATION 


EXECC SOL 2,0 

ALTER 

SUBSCC OPTIONS = K,M,P 

ANY COMBINES 8 REDUCES 
OUTPUT = SUMMARY OF PSEUDOSTRUCTURE 
CONNECTIVITIES 

CASFCC DNA 

BULK DNA 


Table VI. 



IMPLEMENTATION OF A TRAPEZOIDAL RING ELEMENT 


IN NASTRAN FOR ELASTIC-PLASTIC ANALYSIS 

P. C. T. Chen and G. P. O'Hara 
U. So Army Armament Research and Development Command 
Benet Weapons Laboratory, LCWSL 
Watervliet Arsenal, Watervliet, NY 12189 


SUMMARY 


The explicit expressions for an elastic-plastic trapezoidal ring element are 
presented and implemented in NASTRAN computer program. The material is assumed 
to obey the von Mises' yield criterion, isotropic hardening rule and the Prandtl- 
Reuss flow relations. For the purpose of demonstration, two elastic-plastic 
problems are solved and compared with previous results. The first is a plane- 
strain tube under uniform internal pressure and the second, a finite-length tube 
loaded over part of its inner surface. A very good agreement has been found in 
both test problems. 


INTRODUCTION 


In recent years finite element method has been widely used for solving 
complex nonlinear problems and many large-scale general purpose computer programs 
have been developed (ref. 1) . The MARC and ANSYS systems have found wide appli- 
cations, yet they are quite expensive. The piecewise linear analysis option of 
the NASTRAN program provides an algorithm for solving nonlinear problems in mater- 
ial plasticity (ref. 2) . The load is applied in increments such that the stiff- 
ness properties can be assumed to be constant over each increment. However, the 
usefulness of this option is quite limited because only a few elements have been 
implemented. These include rod, tube, bar elements for one-dimensional problems 
and plate elements for two-dimensional plane stress problems. This paper 
describes the implementation of a trapezoidal ring element in NASTRAN for solving 
elastic-plastic problems of rotational symmetry. 

The theoretical basis of our implementation follows the approach first 
developed by Swedlow (ref. 3). A unique relationship between the octahedral 
stress and the plastic octahedral strain is assumed to exist. The material is 
assumed to obey the Mises’ yield criterion, isotropic hardening rule and the 
Prandtl-Reuss flow relations. The explicit expressions for axisymmetric plas- 
ticity are derived. The element stiffness matrix and the stress data recovery 
routines for the trapezoidal ring element are developed. Seven new subroutines 
are implemented and added to the NASTRAN code. 
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For the purpose of demonstration, two elastic-plastic test problems are 
solved. The first is an infinitely long tube under uniform internal pressure. 
The NASTRAN results are in excellent agreement with an exact solution based on a 
finite-difference approach (ref. 4). The second problem is a thick-walled 
cylinder of finite length loaded over part of its inner surface. The NASTRAN 
results are compared with those obtained by a two-dimensional code with the use 
of quadrilateral ring elements (ref. 5) . A good agreement between the two 
results has also been achieved. 


CONSTITUTIVE RELATIONS 


Following the development by Swedlow (ref. 1), the constitutive relations 
to be used in our formulation for solving elastic-plastic problems of rotational 
symmetry will be presented here. In the development, a unique relationship 
between the octahedral stress and the plastic octahedral strain is assumed to 
exist and the use of ideally plastic materials is excluded. The total strain 
components (e r , £ 9 , e z and Y r z) are composed of the elastic, recoverable defor- 
mations and the plastic portions (eP, eg, eP, and yP ). The rates of plastic 
flow, (eP, etc.), are independent of a time Z scale anci are simply used for con- 
venience instead of the incremental values. The definitions of the octahedral 
stress T 0 and the octahedral plastic strain rate feP in the case of rotational 
symmetry are: ' 0 


= (1/3) [(a r -a 0 ) 2 + (a 0 -a z ) 2 + (a z -a r ) 2 + ^rz ^^ 2 
d/3) Uel-kl) 2 + (e ^) 2 + (e ^) 2 + (3/2) (yP z ) 2 ] 


1/2 


CD 

( 2 ) 


where (a r , a 0 , a z , x rz ) are the nonvanishing stress components. 

A unique relationship between x 0 and £§ is assumed and there exists a 
function, M t (t 0 ), such that 

2Hj.(t 0 ) = ijil (3) 


The plastic modulus, Mp(x 0 ), can be related to the slope, E-j., of the effective 
stress-strain ( o-e ) curve by 


1 1 1 
3Mj.(t 0 ) " E t " E 

where E is the elastic modulus and 

* • 

E t = a/e 
a = (3//2)t q 
e = /2 £ 0 + a / E . 


(4) 


(5) 
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The material is assumed to obey the Mises yield criterion and the Prandtl-Reuss 
flow rule. The matrix relationship for the plastic flow in the case of rota- 
tional symmetry is 


[ l * 

k 


l Y 


[C] = i 
E 


rz 
1 r 



( a 


r 


^0 

_ „ 

= [C] « 



°z 


H • 
H 

N 


6T 0 2 MrCT 0 ) 


1 

-V 

-V 

0 

-V 

1 

-V 

0 

-V 

-V 

1 

0 

0 

0 

0 

2 

S r 2 

S r S 0 

S r ; 


s 2 

S 0 ; 


2S T 
r rz 

2S q t 
0 rz 

2S T 
z l rz 


SYM 


4t 


rz -i 


where v is the Poisson’s ratio. 


S r “ ^ 2cr r -a 0~ a z^ 3 » 

S 0 = ( 2a 0-V a r) /3 ’ 
^z — (2 cTz-o^-Oq) /3 . 


( 6 ) 


(7) 


( 8 ) 


For strain-hardening materials, M-p (or E-p) ^ 0, we can obtain the inverse of 
[C] numerically and this procedure is chosen in developing NASTRAN program. For 
ideally-plastic materials, M-p = E-p = 0, matrix [C] does not exist and the NASTRAN 
program fails. However, its inverse [C] _1 still exists and the closed form has 
been derived in Ref. 6. In the axisymmetric case, the explicit form is (ref. 7) 

2G 
A 


[D] = -2*L- 
L J l-2v 


1-v 

v 

V 

0 


SYM' 


1-V 

V 

0 


1-V 


2 a-v) J 


V 


S r S 0 

s 2 
b e 

S r S z 

s Q s 

S r T rz 

s T. 


SYM' 


tz 


St t 2 
z rz rz — 1 


(9) 


103 



where 


and 


A = 3t 0 2 (1 + Hj,/G) 
G = j E/(l+v) . 


( 10 ) 


If we want to remove the limitation that the use of ideally-plastic materials 
is excluded, we have to use Eq. (9) instead of Eq. (7), 

TRAPEZOIDAL RING ELEMENT 

The incremental displacement field employed for the trapezoidal ring 
element are 

Au(r,z) = + 3 2 r + 3 5 z + $ 4 rz , 

Aw(r,z) = 3 5 + 3 6 r + 3yZ + 3 8 rz . (11) 

The transformation from grid point coordinates to generalized coordinates is 

(3) = [r 3q ]{Aq> (12) 


where 


{Aq} = [Auj, Awj, Au 2 , Aw 2> Au 3 , Aw 3 , Au 4 , Aw 4 ) , 

{3} ~ [3p 3 2 > 3^» 3 4 , 3g> 3^> 3y> 3g] • 

and the elements of the inverse of the transformation matrix are 

coefficients of the 3*s in Equations (11), evaluated at the corners of t 


(13) 


the 
the 


element. 


The stiffness matrix is formed in the same manner as that for the aniso- 
tropic elastic element. The final form referred to grid point coordinates is 


where 


[K] - [r^licnrpq] , 

[K] = 27r/r[B] T [D] [Bjdzdr . 


(14) 


[D] , the matrix of material coefficients, is defined by Equation (9). The [B] 
matrix is the same as the elastic case, but now it expresses the incremental 
strains in terms of generalized coordinates 


{Ae} = [B]{3> . 


(15) 
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NASTRAN IMPLEMENTATION 


The NASTRAN implementation for an elastic-plastic trapezoidal ring element 
follows the steps outlined in section 6.8, "Adding a structural element," of 
reference 8. Changes were required in the functional modules PLA1, PLA3, and 
PLA4, which included the writing of several new subroutines. These new routines 
could easily be modeled after the existing code for the linear portion of the 
program. There are two major differences in the nonlinear subroutines. First 
is the new code for the calculation of the material stiffness matrix and second 
is that thermal stresses and element force calculations are eliminated in the 
nonlinear code. 

The changes in PLA1 allows this module to identify the new element as a 
member of the piecewise linear element set and properly initialize the non- 
linear Element Summary and Element Connection Property Tables. Three element 
stress recovery subroutines were added to PLA3: PSAPRG, a driver for stress 

data recovery; PSTRR1 and PSTRR2, phase I and II stress recovery routines. 
Element stiffness calculations in PLA4 require four new subroutines: PKAPRG, 

a driver for nonlinear trapezoidal ring elements in PLA4, PKTRR1 and PKTRR2, 
stress recovery routines which generate stresses for the computation of the 
nonlinear material matrix; and PKTRAP, the stiffness matrix generation routine 
for nonlinear trapezoidal ring elements. 

The computer system available for this work is IBM 360 Model 44. In order 
to add the new code into Link 13, it became necessary to add a new branch to the 
overlay trees to contain the new elements. 


NUMERICAL EXAMPLES 


For the purpose of demonstration, two elastic-plastic problems of rota- 
tional symmetry were solved and compared with other results (refs. 4 and 5) . 

The first is an infinitely long tube under uniform internal pressure and the 
second, a thick-walled cylinder of finite length loaded over part of its inner 
surface. The material properties for both problems are the same. The elastic 
constants are E = 30x10 psi, v = 0.3 and the effective stress-strain curve is 
represented by three line segments connecting the four points in the (e-a) plane, 
(e, 0 ) = (0.0, 0.0). (0.005, 150,000 psi), (0.055, 225,000 psi), (0.1, 225,000 psi). 

EXAMPLE 1. Consider an infinitely long tube subjected to uniform internal 
pressure p. The plane strain condition is assumed. The tube of outside radius 
2" and inside radius 1" has been divided into 25 trapezoidal ring elements. 

The numerical results based on the NASTRAN program have been obtained. For this 
problem, a new finite-difference approach (ref. 4) can be used to generate exact 
solution and to assess the accuracy of the NASTRAN code. Some of the results 
for the displacements and stresses are presented graphically in Figures 1 and 2. 
Twenty-five load increments are used in NASTRAN as shown in Figure 1. The 
radial displacements at the inside as well as outside surface are shown as 
functions of internal pressure. Figure 2 shows the distributions of radial. 
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tangential and axial stress components in a pressurized tube when half of the 
tube is plastic. The pressure required to achieve this state is 0.7378 o 0 based 
on NASTRAN code and 0.7356 o 0 based on the finite-difference solution (ref. 4). 
Both codes indicate that the maximum tangential stress occurs at the elastic- 
plastic interface. As demonstrated in Figures 1 and 2, the NASTRAN results are 
in excellent agreement with those based on the finite-difference approach 
(ref. 4). 

EXAMPLE 2. Consider a two-dimensional elastic-plastic thick-walled cylin- 
der problem as shown in Figure 3. The tube with inner radius (1"), outer 
radius (2") and length (4") is loaded uniformly over a middle portion (2") of 
the inner surface. The mesh generation and the loading for the half of the 
undeformed structure is shown in the figure. This problem was solved in (ref. 5) 
based on a scale loading approach. The first load factor is the upper limit of 
the elastic solution and ten additional increments were needed until one of the 
outside elements becomes yielded. The same load factors were used to obtain 
the NASTRAN solution. Both programs indicate that the sequence in which the 
elements become plastic is 1,5,9,2,13,6,10,3,7,14,11,17,4. Some steps will 
cause more than one. element to become plastic and those elements with effective 
stress a >_ 0.99 a 0 have been considered as yielded. The numerical results for 
the radial displacement at the inside, u a (point 1) and outside, u^ (point 5) 
as functions of internal pressure are shown in Figure 4. The stress components 
at the centroid of one inside element (No. 1) are shown in Figure 5. The effect 
of loading history on the displacements and stresses can be seen from Figures 4 
and 5. A comparison of the results between NASTRAN program and reference 5 
indicates that very good agreement has been achieved. 


CONCLUSION 


An elastic-plastic trapezoidal ring element has been implemented in NASTRAN 
computer program. Its application to elastic-plastic problems of rotational 
symmetry has been demonstrated by solving two thick-walled tube problems. The 
NASTRAN results for both problems are in excellent agreement with the other 
results. 
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NORMAL MODE ANALYSIS OF THE IUS/TDRS PAYLOAD 

IN A 

PAYLOAD CANISTER/TRANSPORTER ENVIRONMENT* 
Karl A. Meyer 

Planning Research Corporation 
John F. Kennedy Space Center, Florida 


SUMMARY 


Special modeling techniques were developed to simulate an accurate 
mathematical model of the Transporter/Canister/Payload (T/C/P) system during 
ground transport of the Inertial Upper Stage/Tracking and Data Relay Satellite 
(IUS/TDRS) payload on the John F. Kennedy Space Center (KSC). The three finite 
element models -- the transporter, the canister, and the IUS/TDRS payload -- 
were merged into one model and used along with the NASTRAN normal mode analysis. 
Deficiencies were found in the NASTRAN program that make a total analysis using 
modal transient response impractical. It was also discovered that inaccuracies 
may exist for NASTRAN rigid body modes on large models when Given's method for 
eigenvalue extraction is employed. The deficiencies as well as recommendations 
for improving the NASTRAN program are discussed. 


INTRODUCTION 


The capability to predict the operational life of a Space Shuttle payload 
is essential to the success of an orbital mission. The amplitudes and approximate 
number of load cycles that will be induced in each individual payload during 
ground transport between facilities must be predictable in advance of each flight 
into space. 

This paper presents a solution to the problem of determining the dynamic 
response of one particular payload: the IUS/TDRS. The TDRS , which will be 
boosted to a geo-synchronous Earth orbit by the IUS, presents a problem in that 
it has a flexible diaphragm tank that contains approximately 1,400 pounds of 
hydrazine fuel. The tank is extremely sensitive to cyclic fatigue. 

Though the original purpose of the study was to determine the dynamic 
response of the IUS/TDRS payload to ground- induced, time-dependent displacements 
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using the NASTRAN computer program, it became apparent that a modal transient 
response analysis (solution 12) was not feasible due to deficiencies in the 
program. 

Instead, it was necessary to perform a normal mode analysis (solution 3). 
Selected results from the NASTRAN program were used as input to a FORTRAN program 
in which time-dependent displacements were imposed and the response was computed. 


THE T/C/P SYSTEM 


The T/C/P system transports Space Shuttle payloads between facilities at 
KSC. It consists of a canister that houses the payload and a transporter on 
which the payload canister is carried. The canister is a steel shell with rigid 
steel frames and two large aluminum doors; it can be used in either a vertical 
or horizontal position. The transporter has a steel flatbed frame that contains 
supporting subsystems such as environmental control, instrumentation and 
communication, fluids and gases, as well as the diesel generators required for 
sel f- propul si on. The frame is mounted on 12 bogie units (6 drive units and 6 
braking units) with 4 tires per unit. 

The IUS/TDRS payload, the subject of this study, is carried with the canister 
in the vertical position. 


NORMAL MODE ANALYSIS 


To perform a normal mode analysis, a NASTRAN model was developed to determine 
the normal mode shapes, eigenvalues, and generalized mass for selected mode 
shapes as well as the total mass, mass moment of inertia, and center-of-gravity 
location of the T/C/P system. It contained 900 grid points in which ASET1 (ref. 
1) bulk data cards were used to reduce the mass Degrees of Freedom to a total 
of 338, using the Guyan reduction technique and the Givens method for eigenvalue 
extraction. A total of 2915 CBAR, CELAS2, CQUAD1, CTRIA1, and C0NM2 elements 
were used in the model . 

Figure 1 shows a NASTRAN undeformed structural plot of the T/C/P system. 
Figure 2 shows a NASTRAN undeformed structural plot of the T/C/P system in which 
the canister was plotted using NASTRAN PLOTEL elements. The second plot 
configuration was used for modal deformed plots to clarify mode shapes. 


Modeling - IUS/TDRS Payload 

The IUS/TDRS payload was modeled using CBAR and C0NM2 elements and connected 
to the payload canister with CELAS2 elements. The IUS/TDRS payload was modeled 
in its own coordinate system using C0RD2R bulk data cards. 
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Modeling - Payload Canister 


The payload canister was modeled using CBAR, CQUAD1, CTRIA1, and C0NM2 
elements. Connections between payload canister and transporter were made using 
CELAS2 elements. The payload canister steel shell was in the basic coordinate 
system, and each canister aluminum door was modeled in its own coordinate system 
using C0RD2R bulk data cards. 


Modeling - Transporter 

The transporter was modeled using CBAR, CTRIA1, and C0NM2 elements. It was 
modeled in its own coordinate system using the C0RD2R bulk data cards. 


DISCUSSION 


Each component of the T/C/P system was modeled in its own coordinate system 
for ease of analysis for future payloads. This modeling technique is particularly 
useful for cases in which the location of the payload relative to the payload 
canister may change and/or the transportation position of the payload canister 
could be either horizontal or vertical. It points out the value of the NASTRAN 
coordinate system bulk data cards: components of a total structural system may 
be translated and/or rotated relative to one another, and only minor changes in 
the NASTRAN bulk data cards are required. 

The NASTRAN PARAM GRDPNT bulk data card was used to locate center of gravity, 
to compute total mass, and to compute mass moment of inertia relative to the 
center of gravity. This output was used as input to the previously mentioned 
FORTRAN program. 

When modeling a plane frame in which the members have open sections and 
substantial depth and are rigidly connected to one another, as in the transporter 
bed frame, the torsional mode eigenvalues are inaccurate unless a specific 
modeling technique is used. This inaccuracy occurs because warping normal and 
warping shear strains (ref. 2) resulting from internal torsional loads are not 
properly accounted for in the model and result in torsional mode eigenvalues 
that are much smaller in value than actually occur. 

To test this fact, two small NASTRAN models of plane frames were analyzed. 
Both models consisted of plane frames made up of I-shaped members; flanges of 
the members were parallel to the plane of the frame. NASTRAN CBAR elements were 
used in both models. The second model had additional members that consisted of 
CBAR elements offset from the neutral axis of the total section to the neutral 
axis of each flange. The PBAR bulk data card for these offset CBAR elements 
consisted only of moment of inertia in the plane of the flange. All other values 
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of this PBAR bulk data card were zero. This modeling technique produces additional 
stiffness terms that increase the torsional rigidity of the structural elements. 
A NASTRAN undeformed structural plot of these models is shown in figure 3. 

A normal mode analysis was performed on both models. The results of the 
first model analysis gave a first torsional mode frequency of 2.716209 Hz. The 
second model first torsional mode frequency was 13.55733 Hz. The error of the 
first model torsional frequency relative to the second model torsional frequency 
was 499%. Figures 4 and 5 show modal deformed plots of the first torsional modes 
of the first and second models, respectively. 

The first model and second model first bending mode frequencies were 44.00233 
Hz and 44.42864 Hz, respectively, which shows that the modeling technique has 
little effect on pure bending mode. It also points out the importance of the 
offset feature on the NASTRAN CBAR bulk data card. Figures 6 and 7 show modal 
deformed plots of the first and second model, respectively. 

This modeling technique was used on the transporter bed frame. Tests were 
run on a similar type frame and results show that this and similar modeling 
techniques produce accurate eigenvalues. 


NASTRAN RESULTS 


The T/C/P system was analyzed using the Givens method for eigenvalue 
extraction in a free body configuration. A Univac 1108 computer was used with 
level 16.0 NASTRAN. 

Table I summarizes the T/C/P system real eigenvalues and frequencies for 
the first six rigid body modes. The eigenvalues and frequencies are not zero, 
thus raising the question of the accuracy of these modes. Also, the force and 
stress output of NASTRAN showed appreciable strain energy at certain locations, 
which raised additional doubts as to the the accuracy of the rigid body modes. 
Because of these questions, the mode shapes were hand calculated at points of 
interest and used as input to the previously mentioned FORTRAN program. 

Table II summarizes the T/C/P system real eigenvalues, frequencies, and 
mode type for the first 12 flexible body modes. No tests related to this T/C/P 
system have been made at this time, so no comment can be made on the accuracy 
of these modes. Figures 8 and 9 show typical modal deformed plots of the T/C/P 
system. 

It should also be noted that the Forward - Backward Substitution time in 
module SMP1 was approximately 23,000 CPU-seconds on a Univac 1108 computer running 
level 16.0 NASTRAN. 
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DEFICIENCIES IN NASTRAN 


It would have been desirable to run a full analysis of the T/C/P system 
(solution 12, Modal Transient Response); however, certain deficiencies became 
apparent as the analysis progressed. These include: 

1. Questionable accuracy of rigid body modes when using Givens method for 
eigenvalue extraction for free body systems 

2. No provisions for non-zero initial conditions relative to modal 
transient response rigid format 

3. No provisions for time-dependent or frequency-dependent displacements 

4. No provisions for checking accuracy of numerical integration for each 
time or frequency step 

Alteration of the NASTRAN program and the use of special modeling techniques 
were considered; however, the time frame of the analysis and the magnitude of 
these deficiencies made these considerations infeasible. 


NASTRAN PROGRAM RECOMMENDATIONS 


The following recommendations are made for the improvement of the NASTRAN 
program for dynamic modal response problems: 

1. Develop a timing equation check of the Forward - Backward Substitution 
method used by NASTRAN. 

2. Check the accuracy of rigid body modes for large problems for all types 
of eigenvalue extraction methods. 

3. Make provisions for non-zero initial conditions as NASTRAN bulk data 
input in the solution of modal response type of rigid formats. 

4. Add time and frequency-dependent displacement equations in the form of 
NASTRAN bulk data cards. 

5. Develop a numerical integration technique (variable step integration) 
that would check each integration step for accuracy. Some allowable 
maximum error should be used as input by the user in the form of a 
NASTRAN bulk data parameter card. 


117 



CONCLUSIONS 


NASTRAN is an excellent tool when used for the extraction of eigenvalues 
and eigenvectors. The use of several methods of eigenvalue extraction and the 
existence of useful types of NASTRAN bulk data cards, such as ASET1, coordinate 
system cards, and PARAM GRDPNT bulk data cards, make it an efficient, time-saving 
program. 

The author concludes, however, that considerable care must be exercised in 
the choice of modeling techniques and that there are deficiencies in the NASTRAN 
program that make a total analysis of some types of problems using modal transient 
response impractical. The above-mentioned recommendations would improve the 
NASTRAN program, making it more useful and efficient for solving modal response- 
type problems. 
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Table I. NASTRAN Rigid Body Modes 


MODE 

NO. 

EIGENVALUE 

FREQUENCY 

(HZ) 

1 

-5.583629 

0.376078 

2 

-3.456232 

0.295884 

3 

-0.273054 

0.083165 

4 

1.667725 

0.205533 

5 

4.727485 

0.346047 

6 

5.503941 

0.373385 


Table II. NASTRAN Flexible Body Modes 


MODE 

NO. 

EIGENVALUE 

FREQUENCY 

(HZ) 

MODE 

TYPE 

7 

723.2128 

4.280093 

IUS/TDRS BENDING 

8 

812.9676 

4.537919 

IUS/TDRS BENDING 

9 

1140.882 

5.375767 

T/C/P SYSTEM 
COUPLED 

10 

1478.755 

6.120238 

TRANSPORTER 

BENDING 

11 

2033.65 

7.177253 

TRANSPORTER 

TORSIONAL 

12 

2330.44 

7.683147 

IUS/TDRS BENDING 

13 

2685.299 

8.247388 

IUS/TDRS BENDING 

14 

2982.683 

8.692079 

IUS/TDRS BENDING 

15 

3236.141 

9.053861 

T/C/P COUPLED 

16 

3517.725 

9.439545 

PAYLOAD CANISTER 
BENDING 

17 

3751.85 

9.748613 

T/C/P COUPLED 

18 

3880.339 

9.914138 

T/C/P COUPLED 
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1 US- TORS PAYLOAO ANALYSIS 

CANISTER - TRANSPORTER - PAYLOAD , CONFIGURATION 3 - RUN I 
KARL MEYER - PRC 1851 - KENNEDY SPACE CENTER 
UNDEFORMED SHAPE 


Figure 1. NASTRAN Undeformed Structural Plot of the T/C/P System 
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IUS-TDRS PAYLOAD ANALYSIS 

CANISTER - TRANSPORTER - PAYLOAD , CONFIGURATION 3 - RUN I 
KARL MEYER r PRC 1251 - KENNEOY SPACE CENTER 
UNDEFORMED SHAPE 


Figure 2. NASTRAN Undeformed Structural Plot of T/C/P System 
with PLOTEL used for Canister Plot 
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PAYLOAD - TRANSPORTER DYNAMIC ANALYSIS 
TEST MODEL 3. WO DO ISC 
ROSELLE HANSON. PRC-1275, KSC 
UNDEEORMED SHAPE 


Figure 3» Test Models 1 and 2 
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PAYLOAD - TRANSPORTER DYNAMIC ANALYSIS 

TEST MODEL », HO 201GC 

ROSELLE HANSON, PRC-1275, KSC 

MODAL OEFOR. SUBCASE l MOOE 7 FREQ, 2.716209 


Figure 4. First Torsional Mode,. Test Model 1 
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PAYLOAD - TRANSPORTER DYNAMIC ANALYSIS 
TEST MODEL 2, NO c018C 
ROSELLE. HANSON, PRC- 1275, K5C 

modal OCrOR. SUBCASE I MODE 7 FREQ. 13.55733 


Figure 5, First Torsional Mode, Test Model 2 
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PAYLOAD - TRANSPORTER OYNAMIC ANALYSIS 

TEST MOOEL i. HO 20I8C 

ROSELLE HANSON. PRC -1 270, KSC 

MODAL DEFOR. SUBCASE 1 MODE 8 FREQ. NH. 00235 


Figure 6. First Bending Mode, Test Model 1 
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IUS-TDRS PAYLOAD ANALYSIS 

CANISTER - TRANSPORTER - PAYLOAO . CONFIGURATION 3 - RUN 1 
KARL MEYER - PRC I SSI - KENNEOY SPACE CENTER 
MOOAL OEFOR. SU8CASE I MODE 7 FREQ. K. 280093 


Figure 8. IUS/TDRS Mode 
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STEADY STATE SOLUTIONS TO 

DYNAMICALLY LOADED PERIODIC STRUCTURES 

Anthony J. Kalinowski 
Naval Underwater Systems Center 


SUMMARY 


The paper treats the general problem of solving for the steady state (time 
domain) dynamic response (i.e., NASTRAN rigid format-8) of a general elastic 
periodic structure subject to a phase difference loading of the type encountered 
in traveling wave propagation problems. Two types of structural configurations 
are considered; in the first type, the structure has a repeating pattern over a 
span that is long enough to be considered, for all practical purposes, as 
infinite; in the second type, the structure has structural rotational symmetry 
in the circumferential direction. Due to the periodic nature of the structure 
and the traveling wave characteristics of the loading, one need only "cut out" 
and subsequently model a typical periodic region of the total structure, 
wherein appropriate periodic boundary conditions (i.e., unknown forces and dis- 
placements are forced equal, except for unknown phase angle, for corresponding 
points on both cuts) are used along the cuts. The paper presents both the 
theory and a corresponding set of DMAP instructions which permits the NASTRAN 
user to automatically alter the rigid format-8 sequence to solve the intended 
class of problems. The new input to a standard version NASTRAN run is a set 
of alter cards, PARAM cards, and direct input matrix (DMI) partitioning arrays 
which are needed for the purpose of partitioning and correspondingly restruc- 
turing the internal NASTRAN mass, damping and stiffness matrices. Final 
results are recovered as with any ordinary rigid format-8 solution, except that 
the results are only printed for the typical periodic segment of the structure. 
A simple demonstration problem having a known exact solution is used to illus- 
trate the implementation of the procedure. 


SYMBOLS 

+ h 

[B] Damping matrix of n xn periodic substructure 

{ F} Total applied force vector 

1 

[I u ] Diagonal unit identity matrix 

I Number of degrees-of-freedom for interior nodes 

[K] Stiffness matrix of n tn periodic substructure 
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L Spatial period of substructure 

r 

L Number of degrees-of- freedom for left cut nodes 

[M] Mass matrix of n*'* 1 periodic substructure 

R Number of degrees-of-freedom for right cut 

t Time 

{11} Displacement vector 

x Spatial coordinate 

to Angular driving frequency 

0, ijj Angle of incident wave 

y Phase constant (complex form) 

y* Phase constant (real part) 

INTRODUCTION 


A periodic structure consists of a number of identical substructures, 
coupled together in identical manners to form the whole system, see for example 
figure la. For such systems, under certain loading conditions, it is often 
possible to treat only one representative substructure in order to obtain the 
general response for the whole system. For example, if the loading is exactly 
the same for all substructures, the latest (and even some earlier) versions of 
NASTRAN can directly solve this class of problem for both static and steadystate 
cases (i.e., rigid formats 1 and 8). In the case of steady state dynamics 
problems (rigid format-8) involving traveling wave propagation type inputs, 
there is a slightly more general loading condition on each periodic substruc- 
ture, namely that the loading on each substructure is identical except for a 
known phase constant y . More specifically, the relation between the applied 
force vector {F} n in the n th substructure and the one, {7}^ in the n^ 1 sub- 
structure is given by 


{F} 


n+1 



(1) 


where y is a known phase constant. For the class of problems addressed in this 
paper, the phase constant is a purely imaginary constant, i.e a . 


y = 0.0 + iy* 


( 2 ) 



and physically refers to the fact that there is no difference in energy loss in 
processing results from one substructure to the next. 

Brillouin (ref. 1) points out that wave motion in periodic systems have 
been studied for nearly 300 years, wherein physicists and electrical engineers 
have worked in this field in problem areas relating to optics, crystals, elec- 
trical transmission liner, etc. (Elachi, ref. 2, provides a comprehensive list 
of 287 references in this field). Applications of this theory to engineering 
structural analysis and solid mechanics type problems is only recent. Refer- 
ences 3, 4, and 5 are typical of analytical solutions to this type problem for 
simple configurations consisting of beams, grillages and plate structures. 
References 6 and 7 represent a significantly more general approach to the prob- 
lem wherein their application of the theory of finite elements enables one to 
solve a much larger class of problems involving rather arbitrary structures 
than one could treat by purely analytical techniques. References 6 and 7 
appear to restrict themselves to the problem of determining the conditions 
(i.e., values of the steady state response freqency w) under which propagating 
or non-propagating free wave motion will occur in the absence of explicit 
external driving forces. 

In the work presented here, periodic structures with explicit external 
driving forces satisfying Equations (1) and (2) are applied to each substructure 
as illustrated in Figure la. 

If the loading and spatial boundary conditions on each substructure are 
the same, except for the phase difference, y*, (i.e„, the loading for two 
typical points in two neighboring substructures separated by the period length 
Lp,are the same except for a multiplying factor of e 1 ’ 1 - 1 *) it follows that the 
response in each substructure is also the same except for the phase difference 
y*. A simple example of such a case is a propagating pressure wave passing 
across an infinitely long ribbed plate as illustrated in figure 2b (the plate is 
in air and no air-structure interaction effects are considered). The propagat- 
ing surface loading wave is given by the formula 


p = p e 1 ( kx + 

r r o 


(3) 


thus the phase difference between any two neighboring substructures is y* = kLp, 
where L^is the spatial period of the periodic system, p Q is the input loading 

pressure amplitude, to is the steady state driving frequency, x is the horizontal 
spatial coordinate, and k is the wave number of the loading wave. 

Other examples of the phase constant relative to a particular example are 
shown in figure 2. In figure 2a, we have a known pressure wave loading propa- 
gating parallel with the axis of the ribbed cylindrical shell; here, the phase 
constant y*is analogous to the figure 2b example and needs no further explana- 
tion. In figure 2c, the incident wave is incident at an oblique angle 0 and 
is incorporated into the formula for the phase constant given in the figure lc. 
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There is a special case where the ends of the structure do not extend to 
infinity (i.e., the ends never meet) but instead are connected cyclically, as 
in figure 3 for example. For such cyclical cases, y* must satisfy an addi- 
tional constraint, namely y* = 2-rr/n where n = 1, 2, 3,* ’ * . For example, in 
the figure 3 case, n = 8, thus y* = tt/4. - 

We limit ourselves to problems having "one-dimensional periodicity", 
whereby this term we imply that only two cuts are needed (we shall refer to 
these as the left and right cuts, see figure lb) to separate the typical sub- 
structure from the system. The response within such a substructure can be, 
however, mul ti -dimensional . The remainder of the paper focuses on the proce- 
dure for obtaining the displacement and stress response within one typical 
block of the periodic system. The typical substructure block can be made up of 
various types of structural elements (including both elements with structural 
damping and nodes with scalar dampers attached) contained within the NASTRAN 
library of elements (e.g., CQDMEM, CQDMEM1, CBAR, CONROD ... etc.). 


SOLUTION FORMULATION 


The solution procedure presented here is very similar to the one in refer- 
ence (6), except for the fact that here we are considering problems with expli- 
cit forcing functions. The first step in the solution procedure is to "cut 
out" the typical substructure from the overall periodic structure as illus- 
trated in figure lb and to subsequently replace the cut nodes with the internal 

forces (OF?) , for the left cut and (F^) for the right cut) that existed at 

AJ II 111 

those nodes before cutting. The displacements at the cut nodes are similarly 
denoted by { U 0 } and {U } where subscripts l and r denote left and right and 

the subscript n denotes the n tn substructure. Since we are only focusing on 

i U 

the results for the n tn substructure, it is convenient to drop the subscript 
n from here on for notational convenience. 

The governing equations of motion for the substructure are first expressed 
in the familiar finite element form 


[M]{U) + [B ] { U) + [K]{U} + (F) 


(4) 


where [M], [B], [K] are the mass, damping and stiffness matrices of the n t * 1 
substructure respectively; {U} is the displacement vector of all nodes of the 
"th 

n tn substructure; {F} is the generalized force fector; (*) = d( )/dt, and bars 
above the variable denote the fact that the variables are complex and that the 

harmonic time response e lwt has not yet been suppressed thus 


{U} = {U}e iwt 


{F} = {F}e iayt 


(5) 
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The next major step is to partition the matrices and vectors of equation 
(4) into left cut unknowns, right cut unknowns and interior unknowns (subscripts 
&, r, i refer to left, right and interior respectively and L, R, I refer to the 
total number of displacement component unknowns for the left, right and interior 
domain respectively; note due to periodicity, L = R. Thus it follows that after 
partitioning we have 


M I M M 

"u 1 “ii £r 


[M] = I £ | M, , | M, r 


H „ | M . I M 
r& 1 n 1 rr 


(L+R+I) x (L+R+I) 




B Jii I B itr 


I- - H 


[B] = 


B U 


B. . B. 
n 1 ir 


rl 


B . B 
ri 1 rr 


[K] 




K *r 

1 ^ 
__J, 

1 ^ 

I K n ■ 
i — i 

K ir 

k h 

1 1 

1 K ri 

K rr 



(1) x (L+R+I) 



( 6 ) 


Note the generalized force vectors {F} has been further decomposed as the sum of 

an unknown force vector, {F 0 }, (which denotes the yet unknown internal forces 

existing at the cuts in the structure) and a known applied force vector, {T 3 }, 
(which denotes all known forces existing within and at the cuts of the periodic 
substructure) . 
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The full periodic structure is cut (left and right cuts), therefore it 
follows that the internal nodal forces normally existing at the cuts now play 
the role of external (as yet unknown) applied forces. 

The special case of an externally applied force appearing at a left or 
right cut requires special attention in that one component of the total force 
vector is due to the externally applied force and the other component is due to 
the internal force at the cut. The external force value on a cut must be 
shared between the generic substructure block being analyzed and its immediate 
neighbor; consequently these end type external force values are divided in half 
(see, for example, the situation in figure lc where node 1 lies on the left 
cut) » 


A further relation that is needed in the formulation relates to the fact 

j_ i- 

that the right end of the n tn generic substructure is the beginning (left end) 
of the n+1 generic substructure, thus from Equation (1) it follows that 


{Fj;} = -e u {F^> 

<v- eU( V 


( 7 ) 


where the minus sign in the first of Equation (7) accounts for the fact that 
internal nodal forces acting as external forces on the right cut of generic 
substructure n are opposite in sign to the internal nodal forces acting on the 
left cut neighboring substructure n+1. 

The next step in the development is to substitute Equations (5) and (6) 

into Equation (4); the subsequent cancellation of e 1a)t permits us to drop the 
bar superscript notation thus arriving at a "reduced form of Equation (4)". 

At this point there are five groups of unknowns, namely (U }, { U 0 } , {U.}, {F^>, 

(F^} 0 The three row partitions of the reduced Equation (4) in conjunction with 
the two Equations (7), provide 3+2=5 corresponding groups of equations to 
balance the five groups of unknowns. Next, we substitute Equations (7) into 
the reduced form of Equation (4), and subsequently employ the third row parti- 
tion of reduced Equation (4) to eliminate the { F?) unknown. Doing these opera- 
tions result in the following set simultaneous equations for the displacement 
unknowns {U^> and {U..}: 
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-O) 2 [MLL] + I -OJ 2 [ML I] + 

iu [BLL] + | ico [BLI] + 

[KLL] | [KLI] 

-u) 2 [MIL] + | -oo 2 [Mil] + 

iw [BIL] + ( iu [B 1 1 ] + 

[KIL] | [KII] 

(L+I) x (L+I) matrix 

where 

[MLL] = [M u ] + Cosh* [M £r J + Cosy* [M^] + [M rr ] 

[BLL] = [B u ] + Cosy* [B^] - wSiny* [M^ r ] + (Siny*/w) [K^] 

+ Cosy* [B r£ ] + [B rr ] + wSiny* [M p£ ] - (Siny*/w) [K r£ ] 

[KLL] = [K u ] + Cosy* [K^] - wSiny* [B £r J 

+ Cosy* [K n ] + [K rr ] + wSiny* [B^] 

[MLI ] = [M £i ] + Cosy* [M H ] 

[BLI] = [B A .] + Cosy* [B^] + o>Siny* [M r - ] - (Siny*/w) [K ri ] 

[KLI] = [K fci ] + Cosy* [K ri ] + wSiny* [B fi ] (9) 

[MIL] = [M. £ ] + Cosy* [M ir ] 

[BIL] = [B u ] + Cosy* [B- r ] - wSiny* [M. r ] + (Siny*/ W ) [K- r ] 

[KIL] = [K u ] + Cosy* [K.^J - wSiny* [B ir ] 

[Mil] = [M ii ] 

[BII] = [B ii ] 

[KII] = [K..] 

At this point, the linear set of complex algebraic Equations (8) can be solved 
for the unknown displacements {U^}, { } . The unknown displacement at the 

right cut, {U r > can be easily computed with the second of equation (7). The 



(1) x (L+I) vector 
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size of the algebraic system is governed by the (L+I) x (L+I) coefficient 
matrix (i.e., matrix multiplying the unknown displacement vector) where L+I 
equals the number of left cut, L, plus interior, I, unknown displacement com- 
ponents. Typically I»L, therefore it is not a very big additional burden on 
the equation solver to include the second of Equation (7) as part of the overall 
system (actually, we add R extra unknowns, (U r >, and R extra equations (where 

R = number of right cut unknowns). Thus in place of Equation (8), we consider 
the slightly larger, but equivalent system of 



-w 2 [MLL] + 
ioj [BLL] + 
[KLL] 

-w 2 (MIL] + 
iw [BIL] + 
[KIL] 

i a) [BRL] + 
[KRL ] 


-to 2 [ML I ] + 
ito [BLI] + 
[KLI ] 

-to 2 [Mil] + 
ito [BII] + 
[KII] 


[ 0 ] 



(L+I+R) x (L+I+R) 
coefficient matrix 



where [BRL] = (siny*/w) [I u ] 
[KRL] = Cosy* [I u ] 
[KRR] = - [I U ] 


( 11 ) 


and [I u ] = diagonal unit identity matrix 

The RI block of the displacement coefficient matrix in Equation (10) above 
is identically zero, thus the bottom R rows of the system of simultaneous equa- 
tions are totally independent of the solution to the top L+I rows. The length 
of the solution vector L+I+R is of exactly the same length of the original sub- 
structure matrix (Equations (4) and (6)), consequently the modification of the 
DMAP instructions becomes simplier because of the fact that one need only 
intercept the logic of the equation solver and replace the existing mass stiff- 
ness and damping matrices with the modified matrices defined by the new coef- 
fient matrix of Equation (10) (and associated new entry definitions from Equa- 
tions (9) and (11)). Since the length of the solution vector is still the same 
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as the original problem before modification, the post processing DMAP opera- 
tions for displacement printout, stress recovery, etc. need not be modified. 

An alternate scheme (although not yet implemented) would be to modify the 
input to the complex equation solver to accept the smaller (L+I) x (L+I) coef- 
fient matrix used in Equation (8) directly. After solving the smaller L+I 
length displacement vector, the full vector (i 0 e», attaching the missing {U^} 

portion) can be formed by expanding it to length L+I+R via the second of Equa- 
tion (7). Finally, stress and displacement results can be processed in the 
usual way with existing DMAP operations. 


RIGID FORMAT -8 DMAP MODIFICATION FOR NASTRAN 


The periodic structure capability described in the previous section can be 
implemented in a standard version of NASTRAN. In particular, the DMAP sequence 
required to perform the necessary operations are listed in Appendix A. This 
DMAP sequence was checked out on an 1108 computer, standard version of level 
15.5 NASTRAN and is introduced in the EXECUTIVE CONTROL deck with the following 
instructions: 

ALTER 138 

(see Appendix A for specific instructions) 

ALTER 139,139 

(replace KDD, BDD, MDD with KDDX, BDDX, MDDX 
within call arguments of FFRD module 
see Appendix A for detailed instruction) 

ALTER 140 

(Conditional print statement, see Appendix- 
A for detailed DMAP instructions) 

These same level 15.5 DMAP instructions can also be applied to level 17.0 
NASTRAN, the only difference being that 

replace ALTER 138 with ALTER 158 

replace ALTER 139,139 with ALTER 159,159 

replace ALTER 140 with ALTER 160 

It is pointed out, however, that the level 17.0 modifications described 
above have not actually been tried although due to the similarity of the change, 
the DMAP sequence is expected to work. 


level 15.5 
implementation 
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It is important to note that we are modifying the standard NASTRAN unknown 
displacement vector coefficient matrix just prior to the entry into the FFRD 
module used for the solution to the simultaneous complex algebraic equations. 
The implication of this statement is that the row numbering scheme for the dis- 
placement vector has already accounted for the fact that single point con- 
straints, multipoint constraints and omitted coordinates have already been 
accounted for. Thus, for example, the length of the {U £ } vector, L, is not 

simply the number of nodes on the left cut times the degrees-of-freedom Der 
node, but rather is less by the amount corresponding to the number of SPC's, 
MPC‘s and OMIT's relating to the nodes along the left cut. Similar comments 
apply to the length of the {U^ } and { U r > vectors. The understanding of the 

above displacement vector length comments must be clearly understood by the 
user before attempting to fill out the input data matrix partitioning vectors 
CV100, CV010, CV001 defined later in this paper. 


INPUT DATA FOR NASTRAN RUN 


The BULK DATA input to a typical periodic structure run consist of two 
basic parts. The first part corresponds to the usual bulk data input cards 
normally required to make a NASTRAN run (e.g., GRID CARDS, ELEMENT CARDS, DAREA 
CARDS, FREQ CARD, DLOAD CARD, etc.); the second part consists of special input 
cards that are explained in the following text. 


PARAM Cards 


These cards are used to enter various matrix coefficients appearing in 
Equations (9) and (11); especially the constants 


Text Variable 

Computer Variable 

Cosy* 

e CMST 

3 

+ 

o 

• 

o 

= Flj^IEG 

-03 2 

= NJ2fMEG2 

-1.0 

e NjZfNE 

-siny*/a) 

e NSMSBj/ 

-(siny*)co 

e nsmst/ 

+1.0 

E ^ note 

s i ny*/w 

e SMSBpf 

(siny*)w 

e SMST0f 


(12) 


0 = zero 
# = letter 
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are read in on standard NASTRAN PARAM cards where y* is the phase angle defined 
in Equation (2) and w is the angular driving frequency in radians per second 
(w = f * 2tt where f = driving frequency in specified on the FREQ card). The 

format for a typical PARAM card is: 


Col's 1-8 PARAM 

Col's 9-16 one of the 8 computer variable names defined by 
Equation (12) 


Col's 17 


24 real part of variable defined in Col's 9-16 


Col's 25 


32 imaginary part of variable defined in Col's 9-16 (only 
non zero entry is for variable FI0MEG) 


Comments 

Strictly speaking, the real part of variable FI0MEG should be 0.0; however, 
for the NUSC Univac 1108, operating with the level 15.5 version of NASTRAN used 
to implement the procedure, an arbitrary small number if entered (say 1.0x10 ) 

in order to avoid a strange system type error message that is printed when 
exactly 0.0 is entered as the real part. The FIJ0MEG variable is only used to 
compute and subsequently print out the internal forces at the cuts after all 
the main calculations for displacement are completed. The mentioned error mes- 
sage probably will not appear if other NASTRAN versions and/or other computer 
systems are used. 


DMI Cards for Matrix Partitioning 


A set of DMI direct matrix input cards are needed to provide the informa- 
tion NASTRAN needs to partition the mass, damping and stiffness matrices. 

Three groups of cards are needed; a column partitioning vector for the left cut 
group of displacement node components, CV100; a column partitioning vector for 
the interior group of displacement node components, CV010; and a column parti- 
tioning vector for the right cut group of displacement node components, CV001. 

A set of row partitioning vectors are automatically generated by the Appendix A 
DMAP instructions. The column partitioning are made up of entries that are 
either 1.0 or 0.0. Since all entries within NASTRAN are assumed to be zero 
unless otherwise specified, the user need only enter 1.0 values in the appro- 
priate slot in each of the above mentioned partitioning vectors. The rules are 
simple and are as follows: 

• Formation of left cut partitioning vector CV100 

Enter a 1.0 in each row number corresponding to each active indepen- 
dent component degrees-of-freedom lying along the left cut. The 
length of the CV100 vector is L+I+R and there should be L 1.0 entries 
(the remaining I+R entries are automatically zero by virtue of not 
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being defined). If the left cut nodal numbering pattern is sequential 
and starts with the lowest node number of the whole system (e.g., node 
1, 2, 3, . . .), then the first L entries of CV100 will be all 1.0 
values. However, if the left cut numbering scheme does not contain 
only the lowest node numbers, but instead the whole system is num- 
bered at random, then the L 1.0 entries will correspondingly be dis- 
tributed throughout the CV100 vector, and the "bookkeeping" involved 
with defining the CV100 vector becomes messy » The user having MPC's, 
SPC s or OMIT's applied to nodes along the left cut must be sure to 
account for these during the process of entering the 1.0 values into 
the partitioning vector CV100. 

Formation of the interior partitioning vector CV010 

Enter a 1.0 in each row number corresponding to each active indepen- 
dent component degree-of-freedom lying on the interior of the struc- 
ture. The length of the CV010 vector is L+I+R and there should be I 
1.0 entries (the remaining L+R entries are automatically zero). If 
the interior nodes are numbered sequentially, (following the same 
sequential pattern used in the CV100 vector), then the middle L+l, 

L+2, . . . L+I entries of the CV010 vector will all be 1.0 values. 
Again remember to account for SPC's, MPC's and OMIT's in the number- 
ing scheme. 

Formation of the right cut partitioning vector CV001 

Enter a 1.0 in each row number corresponding to each active indepen- 
dent component degree-of-freedom lying .on the right cut of the peri- 
odic structure. The length of the CV001 vector is L+I+R and there 
should be R 1.0 entries (the remaining L+I entries are automatically 
zero). If the left cut, interior, and right cut nodes are all num- 
bered sequentially (in the respective order mentioned), then the end- 
ing L+I+l, L+I+2, . . . L+I+R entries of the CV001 vector will all be 
1.0 values. Again remember to account for SPC's, MPC's and OMITS's 
in the numbering scheme. 


DMI Cards Format 


The bulk data cards for the definitions of the partitioning vectors via 
the standard DMI cards is as follows: 

• CV100 vector cards 

-col's .Entry 

first header card 

1-8 DMI 

9-16 CV001 

17-24 0 
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Col 's Entry 

25-32 2 

33-40 1 

41-48 1 

49-56 blank 

57-64 integer value equal to magnitude of (L+R+I) 

65-72 1 

second card 

1-8 DMI 

9-16 CV001 

17-24 1 

25-32 enter integer row number ,say Nl, of first 1.0 entry 

33-40 1.0 for entry Nl 

41-48 1.0 or 0.0 for entry Nl+1 

49-56 1.0 or 0.0 for entry Nl+2 

57-64 1.0 or 0.0 for entry Nl+3 

65-72 1.0 or 0.0 for entry Nl+4 

73-80 Continuation card name, say CONTI, if needed. 

third continuation card (if needed) 

1-8 +0NT1 

9-16 1.0 or 0.0 for entry N1+ 5 

17-24 

25-32 

33-40 

41-48 

49-56 

57-64 

65-72 1.0 or 0.0 for entry Nl+1 2 

73-80 Continuation card name, say C0NT2, if needed. 

fourth continuation card (if needed) 
similar to third continuation card ... etc. 


• CV010 Vector Cards 

These cards are made up analogously to the CV100 vector already described 
above, the only differences being that on the first two cards, replace CV100 
with CV010; also, the number of the first 1.0 entry slot (Nl value in col's 
25-32 of the second CV100 data card) is different. 


• CV001 Vector Cards 

These cards are made up analogously to the CV100 vector already described 
above, the only differences being that on the first two cards, replace CV100 
with CV001 ; also the number of the first 1.0 entry slot (Nl value in col's 
25-32 of the second CV100 data card) is different. 
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Finally, examples of making the partitioning cards is given later in a 
demonstration problem. 


DMI Cards for Merge Operations 

A set of DMI direct matrix input cards are needed to define certain dummy 
matrices which NASTRAN needs to successfully merge certain internal 
matrices within .the Appendix-A DMAP sequence. These input cards will consist 
of a group of eight cards that must always be included for a run. The only 
thing that changes from one NASTRAN run to another is the lengths of these 
dummy arrays. These dummy null arrays are only introduced to avoid DMAP error 
message printouts for Univac 1108, level 15.5 NASTRAN that occurred when the 
lead matrix entry in the standard DMAP MERGE operation is not defined. 


DMI Cards Format 

Col 's Entry 

first header card 
DMI 

LIXLI enter characters, (does not mean multi ply LI times LI) 

2 
1 
1 

blank 

enter integer equal to L plus I 
enter integer equal to L plus I 

second card 

1-8 DMI 

9-16 LIXLI 

17-24 1 

25-32 1 

33-40 0.0 

third card 

same as first card, except 

Col 1 s 9-16 enter LIXL2 

Col's 57-64 enter integer equal to L plus I 

Col's 65-72 enter integer equal to L plus L 

fourth card 

same as second card, except 
Col's 9-16 enter LIXL2 


1-8 

9-16 

17-24 

25-32 

33-40 

41-48 

49-56 

57-64 

65-72 
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fifth card 

same as first card, except 

Col's 9-16 enter L2XLI 

Col's 57-64 enter integer equal to L plus L 

Col's 65-72 enter integer equal to L plus I 

sixth card 

same as second card, except 
Col's 9-16 enter L2XLI 

seventh card 

same as first card, except 

Col's 9-16 enter L2XL2 

Col's 57-64 enter integer equal to L plus L 

Col's 65-72 enter integer equal to L plus L 

eighth card 

same as second card, except 
Col's 9-16 enter L2XL2 


DMI Cards for Unit Matrix Definition 

Finally, a set of DMI direct matrix input cards as needed to define the 
unit matrix [I u ] employed in Equation (11). 


DMI Cards Format 

Col 's Entry 

first header card 

1-8 DMI 

9-16 UMATR 

17-24 0 

25-32 6 

33-40 1 

41-48 2 

49-56 blank 

57-64 enter integer equal to L plus L 

64-72 enter integer equal to L plus L 
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Col's 


Entry 

second card 

DMI 
UMATR 
1 
1 

1.0 

third card 

same as second card, except 

Col ' s 17-24 enter 2 

Col 1 s 25-32 enter 2 

fourth card 

same as second card, except 

Col 's 17-24 enter 3 

Col 's 25-32 enter 3 


etc. 


last (L+l s ^) card 

same as second card, except 

Col's 17-24 enter integer value equal to L 

Col's 25-32 enter integer value equal to L 


DEMONSTRATION PROBLEM 

The use of the DMAP sequence in conjunction with the new PARAM and DMI 
input cards defined in the previous section will perhaps be better understood 
by including a specific example solution that has the features that (1) is a 
small size problem convenient for matrix operation checking and debugging 
purposes; (2) contains most of the main ingredients typical of a representa- 
tive problem, thus the system has mass, stiffness, and damping; (3) the exact 
solution to the problem is known for checking purposes. The problem illus- 
trated in figure 4a meets all of these conditions, and corresponds to a plane 
pressure wave propagating in an infinite acoustic fluid medium. Upon sampling 
the pressure response along any two parallel vertical cuts (separated by the 
horizontal distance Lp) it can be shown that. the response (pressure and dis- 
placement) is different only by the amount e^y* where for this particular 
problem r 


y* = (f)(L p ) Cost (13) 


1-8 

9-16 

17-24 

25-32 

33-40 
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where u> = driving frequency (rad/sec), L p = length (in.) between two parallel 

cuts,.c = compressional wave speed (in/sec); ip = angle of incident wave. The 
solution for the pressure and motion response is sought, for the shaded region 
shown in figure 4a» Only the dotted outlined region is modeled with finite 
elements and is correspondingly shown in figure 4b„ 

In order to exercise this demonstration problem to the fullest extent, we 
will use different types of boundary conditions on all four sides of the fig- 
ure 4b finite element model. 


Boundary Conditions 
Left and right vertical cuts 

Since this problem falls within the class of problems solvable by 
the phase difference type boundary condition, boundary conditions 
specified by Equations (7) are enforced. Here the left and right 
cut forces and displacements are taken as unknowns., 

Top face 

The boundary condition for the top face is different from the left 
and right face in that here we explicitly apply the free field pres- 
sure (converted into equivalent nodal forces). The formula for the 
freely propagating pressure wave is given by the expression 

P(x,y,t) = p Q e l(kr) e la)t ( 14 ) 

P(x,y) 

1 

k = to/c = wave number (in. ) 
r = spatial coordinate normal to direction of wave propagation (in 0 ) 
p Q = steady state pressure amplitude (psi) 

i*fT 

p = pressure (psi) (spatial variation) 
and r is related to the x,y coordinates by the relation 

r = x cosip + y simp (15) 


where 


The y-direction force at upper face node 3 is computed by substituting 
p(x,H) from Equation (14) into the expression 


<3 F I> - 


t 


x=L p /4 


p(x,H) dx 


(16) 
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the force at upper face node 4 is computed similarly by 

x=3L p /4 

V?) = /-p( x > H ) dx 

y J 

x=L p / 4 

and finally the upper face node 9 is computed similarly by 


(17) 


( F?) = 

9 r y 

x=3L p /4 

The demonstration problem is evaluated for the following specific input data: 
p Q = 100, psi 
c = 60000 „ in/ sec 


x=L 


/■ 


-p(x,H)dx 


(18) 


L p = 2 o 0 in. 
H = 2.0 in. 


(19) 


f = w/2rr = 3000o Hz 
ip = 45° 

4 

p = mass density = .000096 1 bo- sec/m 

Substituting the above Equation (19) input constant into Equations (16), 
(17) and (18) results in 

( 3 F®) = -49o8972 e 1 28 » 6378 ° 

(.F?) = -99o7944 e 1 38 » 1837 (20) 

1 y 

( g F®) = -49.9742 e 1 47 * 7297 


An important point must be made regarding loading the final force array 
{F} (Equation (4)) for the periodic structure problem. Observation of the 
Equation (10) loading vector of the new periodic structure problem statement, 

reveals that the left cut applied forces, {F®}, are applied, as normal, to 

the corresponding left cut node; also, the interior applied forces, {F^}, are 

applied, as normal, to the corresponding interior nodes; however, the right 
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• cl 

cut applied forces, (F }, are not applied to the right cut nodes, but rather, 

are first multiplied by the complex constant e“ y \ and then applied to the 
corresponding left cut node, For the demonstration problem at hand, substi- 
tuting Equation (19) into Equation (13) implies that y* = 25,45584°, thus 

e 1w *(,Fj) - -49.9742 e 122 * 27387 ° (21 

y 

therefore, in summary, at node number 3 in the y direction apply a net sum 


force = -49,8972 e 1 28, 6378 .49^742 g i 22 ,27387 
and at node number 4 in the y direction apply a net 


( 22 ) 


force = -99,7944 e 1 38,1837 (23) 

and at node number 9 in the y direction apply a net 
force = 0,0 (i.e., do not apply any force), 

• Bottom face 

The boundary condition for the bottom face could have been selected 
similar to the top face (i.e,, we convert the pressure into equiva- 
lent nodal forces). However, instead we use a slightly more compli- 
cated boundary condition that permits us to introduce a [B] matrix 
entry into the problem. More specifically, the pressure and normal 
velocity along the bottom cut can be obtained from the expression 


p( x ,0) = pcV 


(24) 


where V n is the particlevelocity normal to the wavefront propagation direc- 
tion. Therefore the relation between the V velocity component and the y 


direction resisting reaction is given by 


(Fp 


(p(x,o) *AA) = pcAA Simp 

damping 

constants 


(25) 


Thus, the bottom face boundary condition (simulating an approximate wave 
absorbing boundary condition) is achieved by placing viscous dampers along 
the bottom cut (see figure 4b), wherein the damper constants are defined by 


= pcAA Simp 


(26) 


where AA is an appropriate area factor relating the pressure and concentrated 
force (F. ) . When the wave length of the. incident wave is long relative to 

y 
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the mesh size, one can set AA = the surface element length for bottom surface 
nodes off the cuts and set AA = 1/2 the surface element length for the nodes 
lying on the cuts. Thus for the demonstration problem at hand, C d = 8.14587 

for the middle damper and half that amount for the end cut dampers. 


Preparation of Demo NASTRAN Input Data 

• executive control 

The form of the DMAP instructions presented in Appendix A are general and 
are not problem dependent. The only times the user deviates from the presented 
DMAP sequence is (1) when he switches from one level of NASTRAN to another 
wherein the ALTER statement numbers change; or (2) when he wishes to turn off 
or turn on the intermediate matrix printout switch ISW, defined in the third 
DMAP ALTER statement (ISW = +1 prints all intermediate matrix operation steps 
in addition to a printout of the FORVEC vector which lists the left cut, 
interior, and right cut nodal forces; ISW = -1 no intermediate printout). 

• case control 

The standard CASE CONTROL deck is shown in the APPENDIX A, wherein the 
only thing worth noting is the fact that a DLOAD card is used for the purpose 

of superimposing the two vectors { F„ } and e { F } which are both applied to 

A/ i 

the same left cut nodes (in correspondence with the first partition of Equa- 
tion (10))„ 


• bulk data 

The CDAMP2 cards are used to define the damping constants applied at the 
bottom surface nodes. The CQDMEM membrane elements (and corresponding MAT2, 
PQDMEM cards) are used to define the fluid media, employing the displacement 
fluid element approach described in ref . 8 . The collection of DLOAD, DAREA, 
DPHASE, RL0AD1, TABLED1 cards are used to insert the applied forces defined by 
Equations (22) and (23). The GRDSET card is employed to eliminate the non- 
applicable 3, 4, 5, 6 degrees-of-freedom for the 2-D membrane elements. Stand- 
ard GRID cards define node coordinates and the standard FREQ card defines fre- 
quency of f = o)/2tt = 3000 Hz. The special nine PARAM cards defined by Equations 
(12) are evaluated using the data in Equations (19). In the case of the DMI 
cards for the partitioning vectors CV100, CV010, CV001, the first thing to 
establish for the problem is the sizes L, R, I for the individual partitions. 

It is convenient, for bookkeeping purposes, to number the left cut sequen- 
tially; and to also number the interior sequentially (with the lowest interior 
node number appearing next to the highest left cut node number) and finally 
number the right cut sequentially with the lowest right cut node appearing next 
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to the highest internal node number . Thus for the problem at hand^ L = 2 
degrees-of-freedom/node times 3 left cut nodes = 6; I = 2 degrees-of-freedom/ 
node times 3 interior nodes - 6; and R = L = 6 due to the periodicity of the 
structure., The first 6 entries of CV100 are 1.0, the rest of the entries are 
zero; the nh through 12th entries of CV010 are 1.0, the rest being zero; and 
finally the 13th through 18th entries of CV001 are 1.0, all other entries beinq 
zero. The 8 merging DMI cards LIXLI ••• L2XL2 are defined according to the 
instructions given earlier in the paper and need no further comment. Finally 
the unit matrix DMI cards for the UMATR matrix are defined according to the 
instructions given earlier in the paper. 


Results of NASTRAN Demo Problem 

Selected output results for the figure 4b demonstration problem are pre- 

i^flPPPwnT Y^R TRA Tw dl ^ C ^+ pri ' n ^ 0U x (real and ima 9 inar y P ar t type output format) 
in APPENDIX B. The first part of the NASTRAN printout illustrates the net 

JUS 1 , , ve ^ or results of Equations (22,) and (23 )); it is always good to 

include as a checking feature of the input. 

Next the stress output is printed and refers to the stress computed at the 
° f the element. Note that due to the sign convention difference regard - 
ing stress and pressure (opposite in sign, see ref. 8 for details), the user 
must reverse the sign of stress to obtain the pressure, i.e.. 


p = -a 


xx 


- -a 


yy 


(a P noted that the NASTRAN output format headings are in error 

nf U N^?PMi(°Th tln9 4 . bu S ! n NASTRAN that has remained in practically all levels 
of NASTRAN); the output heading should be read as follows (for Parh pipmpni- 

of stress output): Real Pt. (o )V Real Pt. (a J, Real Pt. (o 1, Img.Pt. 

pt * ( ^ X y ° The resul ts shown in APPENDIX B, Table B-l, 


Img. Pt. (a ), Img 


(a xx ) > 


show the NAZRAN results next to the exact solution, after converting the real 
and imaginary parts into amplitude and phase data. 

r iS eval uate f wi th P( x >y) from Equation (14) at distances 

h8t locate a line (drawn parallel to the wave front) which passes through the 

sol Nt inn °n ^ N ? STRAN results a ^ee well Snh the exact 

solution in both phase and amplitude. It is noted that a very sliqht asvmmetrv 

exists between the NASTRAN phase angle results for element 3 and element 4*. Y 

Probably due to the fact that the boundary conditions on the top and 

thp h £nn U fnr UrfaCSS d ° + 0t r6SUlt in exactl y the same applied nodal forces (e.g., 
t?nn t ?Ll i Ce L W r® no J com P uted using a consistent pressure-to-force applica- 

vcs the displacement shape functions). A similar com- 
ment applies to the relation between the bottom surface forces and left cut 


I° r illust r a t'ive purposes, node 5 had the x displacement SPC con- 
str ^H ec ' to zero and suppose the x displacement of node 4 was forced fthrouah 

cas^L^-r^e^I - q S Tl = J! S R 1 - C L m f n 6 ° f n ° de 6 * In thiS P art ^‘ cular9 
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forces. When the program is run with the detailed print switch triggered 
(ISW = +1), the left and right cut forces are computed from the solution dis- 
placement vector and subsequently printed in condensed format (under the FORVEC 
header). Inspection of these results showed that corresponding left Q and right 
cut internal nodal forces differed by the proper amount, p* = 25.455 . 


CONCLUDING REMARKS 

The DMAP sequence presented in APPENDIX A, permits the NASTRAN user to 
solve a class steady state dynamically loaded periodic structures, wherein the 
internal forces and resulting response at the periodic ends of the typical 

i u* 

repeating substructure are related by a known complex phase relation, e 
The DMAP sequence is general and need not be changed from one problem to the 
next. When using the methodology presented here, one should be extremely care- 
ful that the CV100, CV010, CV001 partitioning vectors are prepared properly. 

For large problems, a preprocessor should be written which automatically gen- 
erates these vectors. One should also be careful to allow for SPC, MPC and 
OMIT cards which compact the solution vector and consequently should be taken 
into account in the preparation of the partitioning vectors. The user is 
strongly advised to study the sample demonstration problem presented here 
before undertaking anymore complicated problems. 

At this point, the DMAP sequence has been checked out for a set of small 
size problems, therefore comments regarding computer run time on large problems 
cannot be made at this time. Some of the DMAP operations can be increased in 
efficiency by employing the matrix operation ADD5 DMAP module in place of 
repeatedly employing the matrix ADD DMAP module. The ADD5 routine was not 
used in the 1108, level 15.5 version of NASTRAN due to the fact that the system 
did not consistently successfully add a string of 5 matrices during certain 
checkout phases of the programming. 

Finally, it is advised that before attempting to exercise the DMAP sequence 
presented in APPENDIX A, one should run the demonstration problem as a bench 
mark check in. If a Univac 1108, level 15.5 version is used, it may be necessary 
to increase the "max-files" to 40 in NTRAN$. 
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APPEND IX -A 

(NASTRAN DEMONSTRATION PROBLEM INPUT) 


N A S T R A N EXECUTIVE CONTROL DECK 


id nusc periodic structure sample: 

APP DISP 

SOL 6»0 

T IME 30 

DIAG 2,3»8»14,15 

DIAS 2 2 

ALTER 138 


$ DMPA INSTRUCTIONS FOR PERIOnTC STRUCTURES (CftnEn By A^J-KALINQMSKTJ 

$ ISW=+1 PRINT ... =-l NO PRINT OF ALL KEY MATRICES PoR PERIODIC STR. 

PARAM //C.N>NOP/V,N.ISW=-H $ 

PARAM //C*N»N0P/V»N.IM1=-1 S 

PARAM //C.N»N0P/V.N.PURSW1=-1 $ 

$ DEFINE ROW PARTITIONING VECTORS WITH DUMMY ADD OPERATION 

$ WHERE CV100 CVOin CVOOl ARE RFAO IN OKI DMT BUI K DATA 

ADD CV100 » /RV1G0 $ 

ADD CVOlOi/RVOlO S _ 

ADD CVOOl »/RV001 $ 

cond mar » isw $ 

MATPRN RV100»RV010»RV001» t// $ 

MATPRN CV100»CV0l0»CVQ01»t// S : 

LABEL MAR $ 

$ CONDITIONAL PRINT OF KDDtMDp, BpQ BEFORE PARTITION APPLIED ( 

COND KAL»ISW * 

MATPRN KDD*MDD»BDD».// * 

$ SAVE ORIGINAL KDD MDO BDD MATRICES FOR PROCESSING FORCES AT CUTS 

ADD KDD,/SAVKPD S ' 

ADD MDD./SAVMDD $ 

ADO BDP./SAVBOD $ 

LABEL KALS 

$ PAR IT ION K»M»B MATRICES 

$ PARTITION LL BLOCK 

PARTN KDD»CV100»RVlfl0/» t »DKLL/C»N»1 $ . 

PARTN MDD.CV100.RV100/»»»DMLL/C*N,1 $ 

PARTN BDD»CV100»RV100/»mDBLL/C>N.1 $ 

S PARTITION LI BLOCK 

PARTN KDP.CVOlOtRVlOO/t mDKLI/C»N,1 $ 

PARTN MDD*CV010»RVl00/»r»DMLI/C»N,l S 

PARTN BD0»CV010»RV100/» * >DBLI/C>N»1 $ 

$ PARTITION LR BLOCK 

PARTN KDDyCVOOl »RVlOO/> t »DKLR/C»N»1 S 

PARTN MDDrCVOOl *RV100/» » »DMLR/C»N» 1 S 

PARTN BDD»CV00l>RVl0Q/» t »DBLR/C>N>1 » 

$ PARTITION IL BLOCK 

PARTN BDD»CV100.Rv010/».»DBlL/CtN»l $ 
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PARTN 

KDDrCVlC0*RV0l0/» » 
MDD»CV100,RV010/» * 

»dkil/c,tt, 

»dmil/c*n» 

1 * 
1 % 



$ PARTITION II BLOCK 





partn 

KDD»CV010»RV010/r t 

»DKH/C*N» 

1 $ 



■fUM 

MDD»CV010»RV010/» » 

7dMiT7c7n7 

1 $ 



PARTN 

BDD»CV01Q*RvQ10/» » 

»DBII/C.N» 

1 s 



$ PARTITION IR BLOCK 





PARTN 

KDD»CV00l#RV010/» » 

»DKIR/C»N» 

1 $ 



ParTn 

MDO#CV001.RV010/» » 

f DMIR/C*N» 

1 $ 



PARTN 

BDD»CV001.RV010/* t 

»DBIR/CiN» 

1 $ 



$ PARTITION RL BLOCK 





PARTN 

KDD> CV100# RVOOl/ 

» r »DKRL/C » 

N»1 

* 


— partn 

MOD » CVlOOr RVOOl/ 

* f »DMRL/C» 

N.l 

$ 


partn 

BDD» CV10C» RVOOl/ 

* * rOBRL/C » 

N.l 

S 


$ PARTITION Rl BLOCK 





PARTN 

KCOf CV010»RV001/» t 

*dkri/c,n> 

1 5 



PART N 

MOD »CV010» RVOOl / * * 

»dmri/c»n» 

1 S 



PARTN 

8DD »CV010# RVOOl / » t 

#dbri/c»n# 

1 S 



S PARTITION RR BLOCK 





partn 

KDD» CVOOli RVOOl/ 

» » » DKRR/C* 

N».t 

s 


PARTN 

MDD» CVOOl# RVOOl/ 

> t »DMRR/C» 

N.l 

s 


PARTN 

BDD* CVOOli RVOOl/ 

»»»dbrr/c» 

N.l 

$ 


SCONDITIONAL PRINT OF PARTITIONS 

cono Kalin * isw $ 


MATPRN DKLL»DMUL»DBLL. .// $ 

m atprn dkli»dmli»obli» »// $ 

MATPRN DKLR»DMLR»DBLR» .// * “ 

MATPRN DKIL>PMIL»PBlL> »// S 

MATPRN DKII*DMII,DBlIr ,// $ 

MATPRN DKIR»DMIR»DBIR» »// $ 

MATPRN DKRL » OMRL » DBrL » t // $ 

MA TPRN PKRI»DMRI»DBRln// S 

MATPRN DKRR»UMRR*DBRR» *// * 

LABEL KALIN $ 

$ 

$ FORM PARTITIONS OF ASSEMBLED MATRIX MDOX» BODX* KDDX 
$ ♦ 

$ FORM LL BLOCK 

1 

$ form amll 

ADD DMlL * BMLK/Di AML l/V* Y »P0NE/V» Y » CMST * 

ADD DlAMLL»DMRL/D2AMLL/V»Y»P0NE/V*Y»CMST $ 

ADD U2AMLL»DMrr/AMLl/V»Y»PONE/V»Y,P0NE * 

PURGE KDD*MDD»BDD/PURSW1 $ 

PURGE D1AMLL»D2AMLL/PURSW1 S 
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"5 FORM ABLL 

ADO DBLL» OB LR/DlBLL/V»Yy PONE/V »Y>CMST $ 

ADD D1BLL * DMLR/D2BLL/V * Y » PONE/V » Y * NSMSTO * 

ADD D 2 BLL * D K LR/D3BLL/ V » Y t PONE/ V t Y > SMSBO $ 

ADD D3BLL • DBRR/D4BLL/V > Y * PONE/V * Y , PONE S 

PURGE D1ABLL»D2ABLL.D3ABLL/PUrSW1 <8 

ADD D48LL » DMRL/053LL/V t Y » PONE/V r Y , SMSTO * 

ADD D5BLL t O KRL/ ABLL/ V > Y t PONE/V t Y f NSMSBQ $ 

PURGE D4BLL * D5BLL/PURSW1 $ 

$ FORM AKLt 

ADD DKLL * OKLR/D1KLL/ V » Y » PONE/V » Y r CMST * 

ADD DlKLL»DBLR/D2KLL/V> Y rPONE/V>Y» NSMSTO < 

ADD D2KLL*0KRL/D3KLL/V»Y»P0NE/V»Y,CMST $ 

ADD D3KLL » 0KRR/D4KLL/V » Y » PONE/V* Y . PONF $ 

ADD D4KLL#DBRL/AKLL/V>Y» PONE/V *Y* SMSTO * 

PU RGE D1KLL t Q2KLL • D3KLL t D4KLL/PURSW1 S 

$ FORM LI BLOCK 

$ 

$ FORM AML I 

A DO DMLI>DMRI/AMLl/V»Y>PONE/V,Y»CMST $ 

$ FORM ABLI 

ADD DBLI »D3RI/D1BLI/V»Y»P0NE/V»Y>CMST i 

ADD D1BLI *DMRl/D2BLl/V»Y»P0NE/V*Y» SMSTO $ 

ADD D2BLI > DKRI/ABH/V»Y > P0NE/V.Y>NSMSBQ $ 

PURGE D1BLI»D2BLI/PURSW1 S 

$ FORM AKLI 

ADD DKLI»DKRI/D1KLI/V#Y»P0NE/V»Y»CMST * 

ADD D1KLI»DBRI/AKLI/V*Y»PQNE/V,Y»SMST0 S 

PURGE D1KLI/PURSW1 $ ~~ 

* FORM IL BLOCK 

* 

$ FORM AMIL 

ADD DMIL»DMIR/AMIL/V*Y»P0NE/V,Y*CMST $ 

* EORMABU, 

ADD DBIL» D9IR/D1BIL/V » Y » PONE/V » Y r CMST * 

ADD DIB IL » DMIR/D2BIL/V > Y > PONE/V » Y » NSMSTO $ 

ADD D2BIL»0KIR/ABlL/VrY»P0NE/V.Y»SMSB0 S 

PURGE P1BIL»D2BIL/PURSW1 $ 

$ FORM AKIL 

ADD DK IL t DK IR/D1K IL/V » Y » PONE/V » Y r CMST * 

"ADD blKlL »DBlR/AKlL/V » Y » PONE/V • Y » NSMSTO $ 

PURGE D1KIL/PURSW1 $ 

1 Form ii Block 

$ AMI I SAME as PMH ABU SAME AS DMH AKII SAMf AS DKII 
S FORM RL BLOCK 
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AUD UMATK » /RBKL/V7Y * SM5B0 5 

AOD UMATR*/AKRL/V,Y,CMST S 

"5" FORM RR BLOCK 

ADD UMaTR»/AKRR/V»Y*NONE $ 

1 THE LR * IR *RI BLOCKS ARE NULL 
$ PRIM BLOCKS BEFORE ASSEMBLY 
COND DEB»ISW $ 

MATPRN AMLL»ABLL>AKLL> .// $ 

MATPRN AMLI*ABLI*AKLI»»// $ 

MATPRN AMIL»ABIL»AKIL» »// $ 

MATPRN ABRL»AKRL»AKRR»»// $ 

PRTPARM //c»n*o$ 

LABEL DEB* 

$ * 

$ NEXT ASSEMBLE BLOCKS WITH MERGE 
$ * 

$ MERGE LL BLOCK 

MERGE LlXLl>»»AKLL»CV100»RV100/KDDXLL/CrN.l $ 

MERGE LIXLI» * »AMLL*CV100»RV100/MDDXLL/C*N»1 * 

MERGE LlXLlr t * ABLL* CV100 »RV10p/BDDXLL/C »N* 1 * 

$ MERGE LI BLOCK 

MERGE LIXL2» r * AKLI »CV010»RV100/KDDXLI/C»N* 1 $ 

MERGE LIXL2»»»AMLI*CV010»RV100/MDOXLI/C»N»1 s 
MERGE LlXL2»»»ABLI»CV010»RV100/BDDXLI/C»Nil * 

1 FORM PARTIAL sum 
ADD KDDXLL»KDDXLI/SUMK1 $ 

ADD MDOXLLrMUDXLI/SLJMMl S 

ADD BDPXLL>BDDXLI/SUMB1 $ 

* PURGE COMPONENTS OF PARTIAL SUM NOT NEEDED ANY MORE 
PURGE KDDXLLfKDDXLI.MDDXLLyMDpXLI.BDDXLL.BDDXLI/PURSwi $ 
$ MERGE IL BLOCK 

MERGE L2XLI r t *AKIL»CV100>RV010/KDDXIL/C»N»1 $ 

MERGE L2XLI * » i AMIL»CV100 »RV0l0/MDDXIL/C »N» 1 * 

MERGE L2XLI t t * ABlL*CV10Q>RV010/8DDXIL/r*N»l $ 

$ CONTINUE PARTIAL SUM 

ADD SUMK1.KPPXIL/SUMK2 $ 

AOD SUMM1»MDDXIL/SUMM2 S 

ADD SUMB1.BDDXIL/SUMB2 $ 

COND JOHN* ISW $ 

MATPRN SUMK1 * KDDXIL* SUMK2 * SUMM1 * MDDXIL//* 

MATPRN SUMM2 * SUMB1 , BDDXlL * SUM02 1 / /% 

LABEL JOHN $ 

1 purge components of Partial sum Not needed any more 

PURGE SUMK1#KDDXIL»SUMM1»MDDXIL»SUMB1.PIDDXIL/PURSW1 $ 

^ MERGE II BLOCK 
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MERGE L2XL2» » >DKII*CV010*RV010/KDDXII/C*N.l * 

MERGE L2XL2> » »DMI I »CV010> RVOlO/MDPXII/C »N» 1 $ 

MERGE L2XL2»»»DBII»CV010»RV010/BODXII/C»N#1 S 

S CONTINUE PARTIAL SUM 

ADD SUMK2*KDDXII/SUMK3 $ 

ADD SUMM2>MDDXII/MDDX $ 

* WHERE MDDX IS THE FINAL MASS MATRIX SUM 

ADD SUMB2.BDPXII/SUMB3 S 

COND BULL.ISW $ 

MATPRN SUMK2tKDnXT I .SUMK3.SUMMj>.MnnXTT//« 

MATPRN MDDX»SUMB2»BDDXII»SUMB3#//$ 

l ab e l BULLS 

S PURGE COMPONENTS OF PARTIAL SUM NOT NEEDED ANY MORE 

PURGE SUMK2 > KQDX 1 1 » SUMM2 1 MDDXI 1 t SUMB2 » BDDXT I/PURSW1 <B 

* MERGE RL BLOCK 

MERGE LlXLI 1 1 t AKRL»rVlQfl»RVQQl /KDnXRt /r »N. i $ 

MERGE LIXLI» > »ABRL#CV100»RV00l/BDDXRL/C»N»l $ 

S CONTINUE PARTIAL SUM 

ADD SUMK3 » KD0XRL/SUMK4 $ 

ADD SUMB3 » BDDXRL/BDPX S 

$ WHERE BDDX IS THE FINAL MATRIX SUM ' 

COND PILL» ISW $ 

MATPRN SUMK3 » KDDXRL . SUMK4 , SUMB3 » BDDXRL//S 

MATPRN BDDX 

LABEL PILLS 

S PURGE COMPONENTS OF PARTIAL SUM NOT NEEDED ANY MORE 

PURGE SUMK3 » KDDXRL » SUMB3 1 BDDXRL/PURSW1 S 
$ MERGE RR BLOCK 

MERGE LIXLI * » » AKRR#CVQ01 » RVOOl/KDDXRR/C* N» 1 S 
$ CONTINUE PARTIAL SUM 

ADD SUMK4 » KDDXRR/KDDX $ ' 

$ WHERE KDDX IS THE FINAL SUMMED MATRIX 
COND BOOK t ISW - $ 

MATPRN SUMK4 » KDDXRR , KDDX » > //S 

LABEL BOOK $ 

S PURGE COMPONENTS OF PARTlAI StJM NOT NEEDED ANY MORE 

PURGE SUMK4»KDDXRR/PURSW1 $ 

s ♦ 

"$ *** NOW ALL THE KDDX BDDX AND MDDX MATRICES ARE FORMED 

$ # 

COND JACK t ISW $ " 

MATPRN MDDX t BQDX t KDDX 1 1 // $ 

LABEL JACKS 

ALTER 139,139 

FRRD C ASEXX * USETD , DLT , FRL , GMD , GOD , KDDX . BDDX » MDDX * » DIt/UDVf , PSF , PDF » PPF/ 
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: . N , DISP/C , N .0 IRECT/V » N t LUSETD/V7N7MPCF1/V#N,SINGLE/V #N*0MT. , 
V » N » NQNCUP/V $ N » FRQSET/C t Y . DEC0M0PT=1 * 


ALTER 140 
COND ABE» ISW $ 


$ PRINT THE FORCE VECTORS USED BY THE EQUIVALENT PERIODIC STRUCTURE 
MATPRN PDF » PSF » PPF » » // $ 

"$ PROCESS THE LEFT AND RIGHT PERIODIC CUT FORCES 


ADD SAVMDC * SAVKDD/DSAV/V » Y » N0MEG2/V » Y » PONE * 


ADD DSAV*SAV0DD/KNET/V» Y f PONE/V» Y • FIOMFG * 
PURGE DSAV/PURSW1 $ 


MPYAD KNET »UDVF» /FORVEC/C »N» 0/C »N»l/c»N»n/C*N*l « 
MATPRN KNET »UDVF»FORVEC 1 1 // $ 

LABEL ABE S 
ENDALTER 


CEND 



DISPLACEMENT s ALL 







CQDMEM 

CQDMEM 


CQDMEM 

DAREA 


6 7 

49.8972 



DMI 

DPHASE 


GRDSET 




UMATR 

23 


6 

3 


4 

3 


3000. 


1 . 

28.6378 


38.18377 

22.27387 



GRID 8 
GRID 


MAT2 10 
PARAM 


PARAM FIOMEG 


345600. 345^00. .0 
.9029168.0 


1.0-20 18849.55 


345600. 
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APPENDIX B (NASTRAN DEMONSTRATION PROBLEM OUTPUT) 

-KtuUENCY = 

3.000000+03 

■ C 0 M P L E X IPAD VErToR 




(REAL/IMAGINARY) 


3 0INT xD. 

TYPE T1 

T2 T3 R1 Rp 

R3 

3 

G * 0 

..t.Q _ „ 

-9.003838+01 .0 .0 .0 

. -4,285621+01 .0 .0 .0 

.0 

.0 

4 

G .0 

-7.844159+01 .0 .0 .0 

.0 

FREQUENCY s' 

c 0 

.0 

3.OO0OOO+O3 

M P L F X ST RE 

-6.169147+01 .0 .0 .0 

.0 

S S E S I N Q U A D R ILATERAL MEMBRANES 

( C 0 D M E M) 


05 


element 

ID. 

1 

2 

3 

4 


real 


xx 


- ST«ESSEs~PJ ELEMENT coordinate system - 
real - • 


. . yy 

■9,069812+01 

ixy 

.0 

lYvr” 

/. 

— y ^ xy 

-4.316528+01 

■9.808248+01 



/. 

.-2^207309+01 

7.895863+01 

• 0 

/. 

-6.209007+01 

9.080797+01 

«_0 

/. 

-4.314245+01 


1ma 9 °yy ~^HEttft-XYircag a 
“4. 310528+01/ » .0 

-2.207309+ 01/. *0 

-6.209007+01/, .0 

-4,3l4g43+01/. .0 


TABLE 6-1 

NASTRAN-EXACT SOLUTION COMPARISON 


elemervE^^^ 

NASTRAN 
|P| .Psi 

EXACT 
|P|» Psi 

NASTRAN 
Phase of p 

EXACT 

Phase of p 

1 

100.53 

100.00 

12.683° 

12.728° 

2 

100.44 

100.00 

25.451° 

25.456° 

3 

100.44 

100.00 

38.180° 

38.183° 

4 

100.53 

100.00 

25.412° 

25.456° 
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Fluid Media 


Plane Pressure Wave 




r=x Cos<Jj+y Simj; 


a) Plane Pressure Wave 
Through Homogenous 
Medi um 



Typical unknown / 
left cut force i- 


boundary 

absorbers 


b) Finite Element Model 
for Dotted Region of 
Fig. 4a 


Fig. 4 Demonstration Problem 


APPLICATIONS OF NASTRAN IN GUST RESPONSE ANALYSIS 

AT NORTHROP 
Ashok K. Singh 
SUMMARY 

A comprehensive gust response analysis .has been performed on a complete 
model of an airplane using the NASTRAN aeroelastic package by the Advanced 
Structural Computer Methods (ASCM) group at Northrop. Earlier the same model 
was used to perform subsonic flutter analysis of the airplane using the 
computer program. Both the random and discrete gust response analyses have 
been performed including the control system dynamics in the problem. On a 
large aircraft gust response analysis including the flexible modes of the 
vehicles is a major design task. On a light weight fighter aircraft the 
analysis is primarily performed in order to study the ride quality and to 
provide the frequency of exceedance curves for the control surface hinge 
moments and some selected dynamic loads for static and fatigue analysis. 

INTRODUCTION 

The NASTRAN finite element program has been used at Northrop since 1972. 

The ASCM group has been actively evaluating and exercising the various NASTRAN 
dynamics analysis features and the aeroelastic package for several months. 
Integrated NASTRAN structural analysis combining static and dynamic analysis, 
e.g., flutter and gust, is the planned goal at Northrop. In order to achieve 
this end, a common structural model throughout an engineering project must be 
used. This practice is also expected to minimize the use of inconsistent struc- 
tural data and unnecessary data handling among the various engineering disci- 
plines . 

Symmetric response analysis is evaluated in two applications to the complete 
airplane aeroelastic model. The analyses are random response to atmospheric 
turbulence and transient response to a discrete gust. The structural and aero- 
dynamic models are the same as that used in the flutter analysis as presented 
in Reference 1. 

The gust response analysis by NASTRAN can only be performed with the aero- 
dynamic forces computed by the Doublet-Lattice Method. Supersonic gust response 
capability has not been provided. Spacewise variation of the gust velocity is 
not allowed in the present NASTRAN formulation but should be considered for 
future development. The gust velocity normal to the free stream velocity is 
taken as an additional source of downwash in the computation of the aerodynamic 
forces. The standard forms of power spectrum of the atmospheric turbulence 
available in NASTRAN are Von Karman and Dryden. However, provision is made to 
use any other form of spectrum by the means of tabular input. In the case of 
the random gust analysis, the frequency response, power spectral density, root 
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mean square value and the frequency of zero crossing, N , of the response are 
output. A post processor can be easily written in order to compute the A and 
frequency of exceedance N(y) of the response. A procedure of weighting the 
aerodynamic forces to match measured data on each panel must be developed. 

The discrete gust analysis is performed by the Fourier method. First, the 
time varying loads are transformed into the frequency domain by Fourier series 
or Fourier integrals. Second, the responses computed in the frequency domain 
are converted into the time domain by inverse Fourier transform methods. 

Three approximate methods are available in order to evaluate the inverse. 

A flight control system is incorporated in the NASTRAN model that utilizes 
a Ride Improvement Mode System (RIMS) and pitch Control Augmentation System 
(CAS) . A dedicated accelerometer and pitch rate gyro near the pilot station 
are chosen to measure acceleration and pitch rate which are fedback as the 
control signal in order to actuate the f laperon and stabilator (Figure 1) . 

To obtain satisfactory ride qualities during low altitude high speed 
flight, extensive effort has been expended in the development of a ride improve- 
ment system. The system tries to maintain a constant value of lift for changes 
in angle-of- attack due to turbulence. This is performed by sensing the load 
factor at the pilot station and using that as a control signal to command the 
high rate flaperons, so that the flaperons can minimize the turbulence-induced 
incremental load factor at the pilot station. 

Operating in parallel with the ride mode is the CAS which uses a blend of 
load factor and pitch rate to maintain aircraft stability while trying to 
minimize uncommanded pitch rate and load factor. 

With the advent of the control configured vehicle, today's control system 
engineer must have a thorough knowledge of aeroservoelastic behavior of the 
flying machine. The design of the filters in the feedback loop cannot be 
completed without the knowledge of elastic and aeroelastic characteristics of 
the modern aircaraft. A new engineering discipline of aeroservoelasticity 
is emerging which will play a prominent role in the early design phases of an 
integrated control system. 


SYMBOLS 

A Ratio of root-mean-square value of load to root-mean-square 

value of gust velocity 

b n b Intensity parameters in the expression for probability of a_ 

H. (co) Frequency response due to the gust excitation 
ja 

N Average number of zero crossings with positive slope, 

per unit time 
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N(y) Number of exceedances of the indicated value of y per unit time 


P P 
1 , 2 


Fractions of total flight time in non-storm and storm turbulence 
respectively 


S (w) Power spectral density of gust velocity 
S (y>) Power spectral density of response 
o r Root-mean-square value of response quantity 

° a Root-mean-square value of gust velocity 

<*> Circular frequency 

co Cutoff frequency beyond which aeroelastic responses are 

c no longer significant in turbulence 

NASTRAN GUST RESPONSE ANALYSIS 

The gust response analysis is performed on a complete aircraft in the 
following steps: 


o Random response to Von Karman gust spectrum without 
control system interaction 

o Random response to Von Karman gust spectrum with 
control system interaction 

o Transient response to a discrete gust with control system 
interaction 


RANDOM RESPONSE ANALYSIS 

A NASTRAN beam element model of a complete airplane was used to perform 
the gust response analysis. The airplane with a tip store, launcher rail, 
wing, flaperon, fuselage, fin with rudder and horizontal stabilizer was modeled 
as finite beam elements as shown in Figure 2. The store and launcher rail 
assembly was tied to the wing tip by rigid elements which may be modified to 
possess elastic properties. The wing root and fin root flexibilities were 
modeled by lumped springs, which may be made more complex as the finite element 
model of the airplane is developed. The horizontal stabilizer root stiffness 
is a general element accounting for the spindle and the actuator assembly 
flexibilities. Mass properties were input on lumped mass element cards. 

A doublet-lattice finite element program was used to represent the aero- 
dynamics of the vehicle as shown in Figure 3. The wing with the launcher rail, 
fin with the rudder and horizontal stabilizer are represented as lifting sur- 
face elements. The fuselage is represented as slender body and interference 
elements. In the present analysis, the aerodynamic induction effect among all 
the elements is considered. The wing, horizontal stabilizer and fin are 
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divided into 145, 20 and 20 micro lifting surface elements, respectively. The 
fuselage is divided into 14 slender body elements. There are 11 interference 
elements on the fuselage. 

A complicated network of linear spline functions is used to relate the 
modal deflections to each of the aerodynamic element deflections. Five 
distinct splines were used for the wing, rail and flaperon panels, two for the 
horizontal stabilizer, three for the fin with rudder, two for the fuselage 
and one for the store. 

The flight control system is represented in the model as 9 extra points 
and 2 scalar points with their coefficients in the mass, damping and stiffness 
matrices in order to represent the filters in the feedback loop. The two 
scalar points are the relative rotations of the flaperon and the stabilator. 
Constraint forces in the equations of motion due to the control laws are 
introduced by the Lagrange multiplier technique. An accelerometer and a 
pitch rate gyro are located near the pilot station in order to measure the 
airplane responses. The measured data are fedback in order to activate RIMS 
and CAS laws, which control the aircraft response. 

The unsteady aerodynamics for Mach 0,8 and sea, level is generated for eight 
reduced frequencies using the Doublet-Lattice Method available in NASTRAN. 
Symmetric flight condition is considered in the analysis. Three symmetric 
rigid body modes and twenty-two elastic modes are selected to generalize the 
aerodynamic forces. Only the modal method of aeroelastic response computation 
is available in the program. Generalized aerodynamic forces at other inter- 
mediate reduced frequencies are computed by means of a linear spline interpo- 
lation routine. 

The gust response analysis of a light weight fighter considered in this 
problem is primarily performed in order to study the ride quality of the vehicle. 
In order to increase the survivability of modern fighter aircraft, new emphasis 
is being given to the capability of low-altitude high-speed penetration. For 
such an aircraft, the low-wing loading /high-lift curve slope which maximizes 
turn rate capability and maneuverability essential for survival in air combat 
also tends to deteriorate the ride quality during high-speed penetration . This 
can lead to reduced mission success in attacking heavily defended ground 
targets, or in the worst case even mission failure. 

For the evaluation of ride quality through turbulence, A is used. A is 
the root-mean-square (rms) of the response divided by the rms gust level in 
feet per second (fps) as defined in Equation 1. 


A = 


f Wc H. (w) 2 S (cu)d<« 
Jo ja a 


0 

r 



S (u))d^ 
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CT 
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( 1 ) 
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To simulate turbulence the Gaussian Von Karman model was used with an rms 
gust intensity level of one fps (References 2 & 3), The scale of turbulence 

used in the analysis is 500 feet. In order to compute the A of the pilot 
acceleration response a cutoff frequency of 30 Hz was used. 

The characteristic frequency, N , is the radius of gyration of the 
response power-spectral density curve with respect to zero frequency (Equa- 
tion 2) , 


“(u>/2tt)®S (u))d«> 

(rvH 


NASTRAN computes a and N by the solution of the airplane equations of 
motion. A post processor may be written in order to compute A and load 
exceedances N(y). Frequency of exceedance, N(y), is the number of exceedances 
of y per unit time or distance flown, where y is any response quantity (Equa- 
tion 3) . 




N(y) = N 

o 


The first fuselage bending mode, whose frequency is 9.5 Hz, is shown in 
Figure 4. Note that the pilot station coincides with the forward node point. 
Transfer functions and power spectral densities of the aircraft with and with- 
out the active controls are given in Figures 5 thru 12. Most of the response 
at the pilot station without the active controls is due to the short period 
and the first wing bending modes. When the RIMS and CAS are incorporated in 
the equations of motion, the response due to the short period mode is markedly 
lowered__but some of the high frequency responses are amplified. The net result 
in the A of the pilot station acceleration is a 37% reduction due to the control 
system dynamics, while the N values are higher. If the visceral response of 
the pilot is only frequency Sependent, the ride quality is significantly 
improved by the RIMS and CAS system used in the analysis. Most of the improve- 
ment in the ride quality is due to the RIMS interaction alone. 

The restart capability of NASTRAN has been used to study the ride quality 
by varying the gains in the feedback loops. Similar restarts may be made to 
vary the scale of turbulence or the gust spectrum at a fraction of the cost 
of the parent analysis. 


exp 
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TRANSIENT RESPONSE ANALYSIS 

Due to the poor ride quality of the modern fighter aircraft a great 
interest has been generated in time domain analysis in parallel with the 
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frequency domain analysis in order to evaluate the handling quality. Transient 
analysis by a Fourier method is available in NASTRAN. The loads defined as a 
function of time are transformed into the frequency domain by Fourier transform 
methods . 

In the present analysis a (1-Cosine) gust profile with a critical gradient 
of 12.5 chord is used. The maximum gust velocity in the profile is 50 ft/sec. 
The NASTRAN restart capability was used to tune the gust response by varying 
the critical gradient. 

Using the Fourier method, the single gust profile is replaced by a series 
of pulses with a period of 20 seconds. The forcing function is zero for some 
time interval to allow for the decay of the responses. In order to evaluate 
the inverse transform equal frequency intervals, method 0, is used. 

The time histories of the relative rotations of the flaperon and the 
stabilator are presented in Figures 13 and 14, respectively . The time history 
of the pilot station displacement is shown in Figure 15. Most of the pilot 
station response is due to the short period mode with a first wing bending 
mode contribution superimposed upon it, as shown by the magnified view in 
Figure 16. The responses are well decayed before the next gust pulse hits 
the aircraft. 


CONCLUSIONS 

This paper shows that NASTRAN is an extremely effective tool for aero- 
elastic analysis. 

A subsonic gust response analysis has been performed in order to evaluate 
its usefulness in the flight vehicle system design. In an earlier evaluation, 
NASTRAN flutter analysis capability has been found to be very satisfactory 
as it provides several state-of-the-art methods and also saves considerable 
amount of man-hours by avoiding duplication of the structural model by the 
static and dynamics group (Ref erence 1). In the pre-NASTRAN era it has been 
a normal practice in many aircraft companies to model the vehicle structure 
separately in the stress, flutter and gust response groups. The aerodynamic 
model is also duplicated within the groups. A consistent and systematic method 
of incorporating control system dynamics in the various engineering disciplines 
is also lacking. With the advent of the control configured vehicle, a new 
engineering discipline of aeroservoelasticity is emerging in a dominant 
engineering role in the early design phases of the aircraft structure and 
integrated control system. NASTRAN aeroelastic package with its integrated 
stress, flutter, gust response, and control systems interaction is a very 
powerful structural analysis tool. It is well tailored for interactive graphics 
environments with data base management systems. With the emphasis on aero- 
elastic tailoring and structural optimization in the aircraft design, NASTRAN is 
amply ready to play a central role. 

The gust response analysis performed in this paper on a light weight 
fighter has given a great insight in the active control system design. The 
use of RIMS and CAS systems reduce the If of the pilot station acceleration 
responses by 37%. The random gust response due to the short period mode is 
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significantly lowered. However, the response due to the structural modes either 
remains unchanged or some of the higher frequency modes are attenuated. This 
means that the present RIMS and CAS feedback loop gains and filters have to be 
modified to suppress the high frequency structural responses. This is an impor- 
tant revelation, which would have gone unnoticed in a customary rigid airplane 
control system analysis. Use of the active control system aggravates the 
dynamic loads on the flaperon, stabilator, etc., structures. 

The transient response analysis of a fighter aircraft is also performed 
by NASTRAN to satisfy a point of view, which wants to look at the time history 
of the pilot response in addition to the frequency related behavior studied 
by the random gust analysis. It is quite possible that the pilot’s discomforts 
are related to the jerks he feels while flying over the discrete impulses of 
the gust rather than a visceral response based on the frequency contents of 
the gust induced excitation. More work has to be done in this area in which 
NASTRAN can be used as a major design tool. 

In all aeroelastic analyses ,the checkpoint and restart capability should 
be used in order to study the influence of changing the scale of turbulence, 
gust spectrum, structural damping, feedback loop gains and f ilters , etc., on the 
responses. Restart procedure is very cost effective and should be further 
improved . 

Generation of aerodynamic influence coefficients and a procedure to weight 
the aerodynamic forces and- moments on each of the panels should be provided 
in order to match the test data. Additional plot capability should be provided 
to plot all the aerodynamic elements, spline fitted modes and aerodynamic 
pressure distributions. For the design of large aircraft , spacewise variation 
of the gust spectrum should be incorporated and coherency, cross-correlation and 
cross-spectral density should be calculated as gust response analysis. Mode 
acceleration method of response computation should be provided in the rigid 
format. A method of computing shear, bending moment and torque at wing stations 
or fuselage stat ions , et c ., and plotting them as response quantities is also needed. 
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Figure 1. Longitudinal Control System Concept 
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Figure 2 


Complete Airplane Structural Model - NASTRAN 
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IMPLEMENTATION OF NASTRAN ON THE 
IBM/370 CMS OPERATING SYSTEM* 
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SUMMARY 


The NASA Structural Analysis (NASTRAN) computer program is 
operational on the IBM 360/370 series computers. While execution 
of NASTRAN has been described (ref. 1) and implemented under the 
virtual storage (VS) operating systems of the IBM 370 models, the 
IBM 370/168 computer can also operate in a time-sharing mode under 
the virtual machine (VM) operating system using the Conversational 
Monitor System (CMS) subset. This report describes the changes re- 
quired to make NASTRAN operational under the CMS operating system. 


"The views and conclusions contained in this 
document are those of the contractor and should 
not be interpreted as necessarily representing 
the official policies, either expressed or implied, 
of the United States Government." 
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INTRODUCTION AND BACKGROUND 


M.I.T. Lincoln Laboratory first obtained NASTRAN in April 1974, when it 
purchased Level 15.5 from COSMIC for use on its newly installed IBM 370/168 
computer. Minor modifications were made in order to make NASTRAN operational 
under the virtual storage concept of the IBM 370 system. However, much of 
the source code did not match the special executable load module code that 
had been received for direct loading onto IBM 3330 disk drives. 

The Laboratory's NASTRAN capabilities were upgraded to Level 15.5.3+ in 
November 1975 with the assistance of NASA Goddard Spaceflight Center. The 
new source code and new executable load module code was in complete agree- 
ment. This Level 15.5.3+ of NASTRAN has adequately served the computational 
needs of the Laboratory since that original installation. 

However, the Laboratory's computer configuration and operational schedule 
restrict the original version of NASTRAN to nighttime operation under the 
virtual storage (VS) batch-processing environment. The daytime operating 
mode (when engineering users are present) consists of an interactive time- 
sharing environment of up to 170 users under the virtual machine (VM) opera- 
ting system using the Conversational Monitor System (CMS) subset. 

It was desired to have NASTRAN available for general use in the inter- 
active (CMS) as well as the batch (VS) environment, with commonality of input 
files. From an engineering user standpoint, the availability of NASTRAN in 
both environments is highly desirable. The CMS environment can be used for 
input data syntax checking, plotting of input and/or output data, and execu- 
tion of relatively small analyses; while the VS environment can be used for 
execution of large analyses that require substantial computer time. The CMS 
environment is also useful to the programmer for maintenance and development 
of the many NASTRAN subroutines and functions. 

For these reasons, the decision was made to develop a CMS/370 version of 
NASTRAN, working from the Laboratory's Level 15.5.3+ source. Compatabil ity 
between the source files, the source listings, and the executable code under 
both the VS and CMS operating systems was required. Our primary objective 
was to make NASTRAN operational with minimal changes in computer coding. 

This report describes the changes which were made to the NASTRAN computer 
program to develop a time-sharing version under 370-CMS and a compatible 370-VS1 
batch-processing version. In this report, the following terminology is used: 

1. CMS- NASTRAN: the version of NASTRAN developed at M. I .T. 

Lincoln Laboratory to operate under VM/370 (Virtual Machine 
facility for the IBM S/370) using the CMS (Conversational 
Monitor System) subset of VM/370. 
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2. VSl -NASTRAN: the version of NASTRAN developed at M.I.T. 

Lincoln Laboratory to operate under the VS1 (Virtual 

Storage) operating system for the IBM S/370. 

The following definitions are based on the hardware and software devel- 
oped by IBM for the implementation of the virtual storage concept. 

Address Space - the set of memory addresses used by a 
program. VS1 allows a maximum total address space 
of 16 megabytes. CMS allows each user a possible 
maximum address space of up to 16 megabytes (limited 
at MIT-LL to 2 megabytes). 

Page - a subdivision of the address space , 2K bytes on 
VS1 and 4K bytes on CMS. 

Real Memory - the set of memory addresses which are 
physically available on the CPU. (3 megabytes 
on the Lincoln computer). 

Virtual Memory - the address space which can be 
addressed by a relocating CPU. Physically, 
the virtual memory exists on a direct access 
storage device, and although a program may 
reference virtual memory in a random fashion, 
the information must be transmitted from virtual 
memory (direct access) to real memory one page 
at a time . 

CMS DIFFERENCES AFFECTING NASTRAN 

The CMS subset of the VM operating system is noticeably different from 
its VS operating system counterpart, in several areas: 

(1) Input/output file structure, 

(2) Load module formation, and 

(3) Memory management and open core concepts. 

These differences between CMS and VS required design decisions regarding 
implementation of NASTRAN under CMS. In all cases, we attempted to introduce 
minimal coding changes into the Level 15.5.3+ version of NASTRAN and tried 
to remain consistent with the original design philosophy regarding NASTRAN 
operation on IBM computers. 

The following sections will deal in specifics concerning the above- 
mentioned areas of difference and will describe our solutions to the problems 
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encountered. The final section will mention some discovered coding errors 
as well as the implementation of a plotting package using the SC4060 plotter. 

Input/Output File Structure 

Execution of NASTRAN under the VS operating systems requires a large 
number of standard Job Control Language (OCL) cards. Modification of the 
Procedure (PROC) card that precedes these basic JCL cards allows the user 
flexibility to define space allocations, physical units, and program libraries 
to be searched. 

The CMS operating system has NO Job Control Language cards. File defini- 
tions are performed by issuing FILEDEF commands immediately prior to execution 
time. The FILEDEF command defines the physical device type and the character- 
istics of the file. All CMS files are always stored and retrieved in 800 
byte blocks. We have chosen 6400 byte blocks for the NASTRAN files - a size 
which will physically fit two blocks per track on an IBM 3330 disk drive as 
well as be wholly compatible with CMS block sizes. These FILEDEF commands 
are stored into an EXEC file which can be invoked by a user along with other 
commands to load and execute CMS-NASTRAN. (refer to Appendix A). 

In order to minimize the coding changes required in GNFIAT (the program 
that generates the file allocation tables), it was necessary to provide a 
list of acceptable file names as input to the program. Under the VS operating 
system, acceptable names are chosen from the JCL cards included in the data 
deck. In CMS, the first card in the IOLIST file indicates (412 format) the 
number of acceptable (permanent, primary, secondary, and tertiary) file names, 
and the names themselves follow (refer to Appendix B). 

Since file chaining is accomplished by CMS, the SPACE limitations needed 
for VS operations are no longer pertinent under CMS. For this reason, all 
IOLIST files were allowed to expand in additional 6400 byte blocks as often 
as needed to provide file space. The VS1-NASTRAN version uses the same file 
extension scheme outlined (ref. 2) for NASTRAN on the IBM 360-370 operating 
system. 

The POINT macro differs in its operation under CMS from that under VS 
and thus coding changes were necessary for successful operation of BSAM 1/0 
processing in CMS-NASTRAN. 

Load Module Formation and Usage 


The most significant difference concerning load module formation and 
usage under the CMS operating system is that no provision exists in CMS for 
an overlay structure. The NASTRAN load module generation under VS attempts 
to minimize the core requirements of each control routine LINKNSii, ii =01, 
..., 14 by use of overlay structures. Program core space in a virtual machine 
no longer becomes so critical, and thus provisions for overlay structure under 
CMS were felt not to be important. 
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As a simple first approach, the load module generation of all fourteen 
control routines LINKNSii was carried out with no overlays. While LINKNS01 
required 660K for storage, the next largest, LINKNS13, required only 530K. 
NASTRAN, the super- link module, needs 45K and the system consumes 129K. Allow- 
ing for a reasonable work space in open core and providing FORTRAN with 
sufficient buffer space, CMS-NASTRAN thus requires a minimal core of 1200K. 

Our present account structure provides us with 2048K which we fully utilize 
while running CMS-NASTRAN. However, we are eyeing the possibilities of 
reducing core requirements to 1024K by splitting the larger LINKNSii control 
routines into multiple modules. 

Load module formation under CMS also differs from the process of linkage 
editing used by VS. The linkage editor under CMS first searches the user's 
directory for TEXT files with the same names as CSECT entry names and then 
searches user-specified libraries called TXTLIBs in a specified order. Unlike 
parti tioned-data sets (PDSs) on VS where member names are user-specified, the 
CMS main entry in a TXTLIB consists of the name of the first CSECT in the 
program, irrespective of the original name of the compiled program. For this 
reason, FORTRAN subroutines in NASTRAN containing only BLOCK DATA statements 
and named COMMONS have the CSECT name of the first named COMMON when compiled 
and thus that entry name when stored in TXTLIBs under CMS, but have the 
assigned member name when stored in OBJECT PDS's under VS. 

In VS, the linkage editor has input cards - such as INCLUDE and INSERT - 
which specify specific members to be included in a module from a PDS and 
specify specific placement of a CSECT within a load module. It was this ability 
of the linkage editor to have the user specify placement within a module which 
enabled the developers of NASTRAN to implement the open core concept on IBM 
batch systems. 

In CMS, there is no linkage editor, per se, but rather a loader with two 
user functions associated with it - LOAD and INCLUDE. The loader builds its 
own CSECT list which can be resolved from only two places: TEXT files or 

TXTLIBs. All TEXT file names and entry names in TXTLIBs must therefore be 
CSECT names. Since the loading process is primarily a search process, the 
order of CSECTs in a module will depend on the order in which programs were 
encountered in the search. This absence of user-specified ordering in CMS 
loading required the removal of certain* CSECTs (CONMSG, MAPFNS, EJDUM2, and 
some named COMMONs) from TXTLIBs and their storage in TEXT files with file 
names different from CSECT names. These CSECTs will be unresolved externs 
after the LOAD function but will then be put at the end of the module in a 
user-specified order by means of the INCLUDE function, making sure EJDUM2 is 
last since it represents open core. 

The CMS linkage editor also differs from its VS counterpart in that it 
resolves COMMON CSECTs at execution time rather than at load time. For this 
reason, a way had to be found to force those FORTRAN subroutines containing 
only BLOCK DATA statements and named COMMONS to be loaded with the other CSECTs. 
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This was accomplished by writing additional FORTRAN subroutines LINKiiC that 
did nothing more than CALL the named COMMONS whose CSECT entries existed in 
the TXTLIBs. This forced inclusion at load time. (Refer to Appendix C.) 

A change was also required to LINKNSii and NASTRAN to reference the proper 
entries in the current FORTRAN function library (e.g. IHOERRM, etc.). 

The lengths of all load modules formed in the CMS system are stored in a 
MASTER file. These lengths can than be used to control the open core require- 
ments described in the next section. Load module generation under CMS reserves 
the first 128K (origin hex 20000) for the operating system and CMS. All load 
modules have the user's lowest address origined to 20000, which for NASTRAN 
execution meant loading the super-link NASTRAN at location 20000. In order 
to allow for future growth and/or modifications to the super-link NASTRAN, the 
origin for all LINKNSii entries was taken to be hex location 2AD00. This 
LINK origin is stored in the MASTER file along with the total load module 
lengths of all fourteen LINKNSii modules (refer to Appendix D). 

Memory Management and Open Core Concepts 

Memory management and open core concepts under the VS operating system 
have been very clearly presented (ref. 1) and will not be reviewed here. In- 
stead, the manner in which the CMS operating system must handle these tasks 
will be explained. 

Under CMS, the super-link NASTRAN load module is initially loaded into 
core at origin hex 20000 and execution begins. (Figure 1 depicts address 
space allocation for CMS-NASTRAN execution.) 

1. Within the NASTRAN module, a user request for 
working storage (GETMAIN) will be issued for all 
of available memory. 

2. The NASTRAN module then releases F0URK (16K words 
under CMS) high address space for operating system 
use via the FREEMAIN macro. 

3. All of the GETMAINed area will be managed by 
NASTRAN rather than by CMS for all executions 
of LINKNSii load modules requested by the 
super-link program NASTRAN. 

4. The load modules LINKNS01 thru LINKNS14 have been 
generated at origin 2AD00 and are. loaded with the 
CMS macro L0ADM0D instruction as they are needed. 
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The above memory management scheme was required since the CMS simulation 
of VS LINK, 6ETMAIN, and FREEMAIN macros did not present a duplicate image of 
core to NASTRAN of that which VS presented. Thus, the macro for loading a 
module into core could not be LINK but had to be LOADMOD. This in turn said 
memory had to be obtained prior to and for the LOADMOD. Since the LOADMOD 
forced CMS to obtain buffer space for itself, we were forced to resort to a 
single GETMAIN with total core management performed by NASTRAN so that frag- 
mentation of memory (into NASTRAN blocks intersperced with system blocks) was 
eliminated. 

Since CMS-NASTRAN has now assumed the responsibility of managing memory, 
the "open core" concept of data management must also be assumed by NASTRAN. 
Under the overlay load module structure of VS, "open core" at the end of an 
overlay tree segment was denoted by a dummy named-COMMON section. This named- 
COMMON could contain an array whose length extended into the open core region 
but yet would not destroy or overwrite code or data in other segments of the 
overlay. Unfortunately, under CMS with no overlay structure, we are faced 
with numerous dummy named-COMMONs , each of which must be located near the 
end of the load module so that it will not overwrite other code and/or data 
in the module. The placement of these named-COMMONs is accomplished by 
INCLUDing newly-written FORTRAN subroutines, LINKiiCC, when doing the load 
module generation (refer to Appexdix E) . 

Minor Differences 


Several subroutines contain DATA statements which are used to initialize 
variables which are subsequently modified. The DATA statements were used 
assuming that the programs containing these statements exist on segments of 
the overlay tree in VS, implying that a fresh copy of the subroutine will be 
loaded (and thus reinitialization of DATA variables will occur) each time the 
segment is needed. Under CMS, the subroutine when it is reentered contains 
the last values of the DATA variables and they are not reinitialized, thus 
causing errors. This problem would also occur in a non-overlay VS version of 
NASTRAN. 

The graphics package acceptable to the CMS environment at M.I.T. Lincoln 
Laboratory utilizes SC4060 meta code to generate the meta code for on-line 
Tektronix terminals. For this reason, a SC4060 plotting package was added 
to the NASTRAN package. Since the CMS environment provides disk file facilities 
for graphics as well as standard I/O, restrictions on graphics files in NASTRAN 
were removed So that disk file definitions for PLT1 and/or PLT2 files are 
permissible. 

Conclusions 


The implementation of CMS-NASTRAN has been completed. While operation 
has been limited to test examples and several small production runs, it 
definitely shows promise as a useful program on future projects in which 
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NASTRAN will be needed. Program enhancements (such as the SC4060 graphics 
package) are easily made while working in the CMS environment and can be 
added to the VS-NASTRAN after successful debugging has been completed. 

A new VS1 -NASTRAN is being generated at the present time. It will embody 
the design philosophy of CMS-NASTRAN and, thus, will not be in overlay struc- 
ture. Compatibility of source (except for several assembly language programs) 
and executable code between the two systems should provide for the reliability 
needed to assure identical results in either mode. 

It is anticipated that CHKPNT saves and restarts should be possible be- 
tween the CMS and VS operating systems. This would allow preprocessing (data 
syntax checking, undeformed structural plots, etc.) to be performed and check 
pointed under CMS; restored, executed, and checkpointed under VS; and finally 
restored and postprocessed (including deformed structural plots) under CMS. 
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Appendix A : Partial Listing of NASTRAN EXEC File 

SCONTROL OFT NOMSG 

-START STYPE ENTER FILENAME OF NASTRAN FILE 
BREAD ARGS 

SIF SINDEX NE 1 SGOTO -START 
SFN = SI 

-CONT STYPE ENTER FILEMODE FOR NASTRAN OUTPUT FILES 
BREAD ARGS 

SIF SINDEX NE 1 &GOTO -CONT 


SFM = & 
FILEDEF 

1 

1 

DISK 

■&T01 

NAST 

&FM 

(BLOCK 

80 

P.ECFM 

F 

LRECL 

80 

FILEDFF 

4 

DISK 

FT 04 

NAST 

&FM 

(BLOCK 

133 

RECFM 

Y 

LRECL 

133 

FILEDEF 

5 

DISK 

SFN 

NASTRAN 








FILEDEF 

6 

DISK 

SFN 

OUTPUT 

SFM 

(BLOCK 

133 

EECFM 

F 

LRECL 

133 

FILEDEF 

7 

DISK 

SFN 

PUNCH 

SFM 

(BLOCK 

80 

RECFM 

F 

LRECL 

80 

FILEDEF 

POOL 

DISK 

POOL 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

NPTP 

DISK 

NPTP 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

OPTP 

DISK 

OPTP 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

PLT2 

DISK 

SFN 

4060 

SFM 

(BLOCK 

800 

RECFM 

F 

LRECL 

800 

FILEDEF 

PRI01 

DISK 

PRI01 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

PRI02 

DISK 

PR 10 2 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

PRI03 

• 

• 

DISK 

PRI03 

NAST 

SFM 

(BLOCK 

6400 

• 

• 

RECFM 

F 



FILEDEF 

• 

• 

PRI30 

DISK 

PR 13 0 

NAST 

SFM 

(BLOCK 

• 

• 

6400 

RECFM 

F 



FILEDEF 

PRI31 

DISK 

PRI31 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

PRI32 

DISK 

PR 13 2 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

SEC01 

DISK 

SEC01 

NAST 

SFM 

(BLOCK 

64 00 

RECFM 

F 



FILEDEF 

SEC02 

DISK 

SEC02 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

S ECO 3 

DISK 

SEC0 3 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



FILEDEF 

TER01 

DISK 

TER01 

NAST 

SFM 

(BLOCK 

6400 

RECFM 

F 



LOADMOD 

NASTRAN 












START 



Appendix B : Listing of NASTRAN IOLIST File 


04320301 

PERM POOL 

PERMNPTP 

PERMOPTP 

PERMPLT2 

PRIMPRIO 1 

PRIMPRI 02 

PRIMPRI03 

PRIMPRI04 

PRIMPRI05 

PRIMPRI06 

PRIHPRI07 

PRIMPRI 08 

PRIMPRI09 

PRIMPRI 1 0 

PRIMPRI 1 1 

PRIMPRI 1 2 

PRIMPRI1 3 

PRIMPRI1 4 

FRIMPRI 1 5 

PRIMPRI 1 6 

PRIMPRI 1 7 

PRIMPRI 1 8 

PRIMPRI 1 9 

PRIMPRI20 

PRIM PRI2 1 

PRIMPRI22 

PRIMPRI23 

PRIMPRI24 

PRIMPRI25 

PRIMPRI26 

PRIMPRI27 

PRIMPRI28 

PRIMPRI29 

PRIMPRI30 

PPIMPRI31 

PRIMPRI32 

SECOSECO 1 

SECOSEC02 

SECOSEC03 

TERTTERO 1 
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Appendix C : Listing of Typical LNKiiC FORTRAN Subroutine 


C 

c 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE LNK01C 

SUBROUTINE BPDAED 
CALL XSFA1 
SUBROUTINE IEXDBD 
CALL IFPDTA 
SOBBOUTINE IFXOBD 
CALL IFPXO 
SUBROUTINE IFX1BD 
CALL IFPX1 
SUBROUTINE IFX2BD 
CALL IFPX2 
SUBROUTINE IFX3BD 
CALL IFPX3 
SUBROUTINE IFX4BE 
CALL IFPX4 
SUBROUTINE IFX5BE 
CALL IFPX5 
SUBROUTINE IFX6BD 
CALL IFPX6 
SUBROUTINE IFX7BE 
CALL IFPX7 
SUBROUTINE IFPABE 
CALL IFP1A 
SUBROUTINE AXICBD 
CALL IFP3BE 
SUBROUTINE XGPIEE 
CALL XGPIC 
SUBROUTINE XMPLBB 
CALL XGPI2 
SUBROUTINE XSORBE 
CALL XSRTBB 
SUBROUTINE UMFZBE 
CALL UMFZZZ 
SUBROUTINE XBSEE 
CALL XLKSPC 

RETURN 

END 
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LINKORIG 
LTNKNSO 1 
LINKNS02 
LINKNS03 
LINKNS04 
LINKNS05 
LINKNS06 
LINKNS07 
LINKNS08 
LINK NS09 
LINKNSIO 
LINKNS1 1 
LINKNS 1 2 
LINKNS 1 3 
LINKNS 1 4 


Appendix D : 


Listing of NASTRAN MASTER File 


02AD00 
ODOOA8 
08C580 
092F50 
072030 
094F00 
078218 
078428 
0976B0 
054950 
073AD8 
08D8D0 
05C688 
OAF 908 
07A390 
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: Listing of Typical LNKiiCC FORTRAN Subroutine 


C DUMMY TO FORCE XGPI1 TO COME JUST BEFORE EJDUM2 (OPEN CORE) 

BLOCK DATA 
COMMON DUM4 (100) 

COMMON /GINOX/ DU M6 (163) 

COMMON /SETUP/ DUM2(7) 

COMMON /IFP3LV/ DUM10(104) 

C BEGINNING OF DUMMY COMMONS THAT MARK START OF OPEN CORE. 

COMMON /IFP1X/ DUM9 (371) 

COMMON /ENDSSS/ DUM7 (2) 

COMMON /I FP XX/ DUMB 
COMMON /IFP3ZZ/ DUM 1 1 
COMMON /IFP4ZZ/ DUM12 
COMMON /IFP5ZZ/ DUM 1 3 
COMMON /XCSABF/ CUM14 
COMMON /ESOFT/ DUM 
COMMON /UMFXXX/ DUM1 
COMMON /ESF A/ DUM5 
COMMON /X GP1 1/ DU M 20 ( 5) 

COMMON /EJDUM2/ DUM21 
END 
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IFEMS, AN INTERACTIVE FINITE ELEMENT 
MODELING SYSTEM USING A CAD/CAM SYSTEM 


Spencer McKellip, Todd Schuman and Spencer Lauer 
Sikorsky Aircraft 
Division of United Technologies 


SUMMARY 


General purpose 3D finite element mesh generators, require 
detailed geometric description of the component to be modeled. 

For complex shapes, this process can be very time consuming. On 
the other hand, most general purpose CAD/CAM systems have the 
capability of defining geometry in an efficient, interactive en- 
vironment. This paper describes one method of coupling a CAD/CAM 
system with a general purpose finite element mesh generator. 

The problem is one of generality in appl ication . . Most inter- 
active mesh generators are either tailor made to specific appli- 
cations or are restricted to classes of geometry suchas bodies of 
revolution, dragged shapes, 2D geometry or digitized input. What 
is needed is a completely general purpose, 3D interactive tool 
that can handle any conceivable geometry without digitizing or any 
other restriction. 


A three dimensional interactive graphic system for defining 
geometry can be an expensive proposition for just finite element 
analysis. However such capabilities have been in existence for 
some time in the form of computer aided design and manufacturing 
systems, commonly known as CAD/CAM systems. These systems are 
logical candidates for front ends to a general purpose, 3D finite 
element mesh generator. 


The 

1) 

2 ) 

3) 


system described in this paper consists of three programs: 

CAIDS, A Sikorsky developed CAD system. 

IFEMS, A Sikorsky developed interactive processor to 
couple a CAD system with a mesh generator. 

PWAMESH , A Pratt & Whitney developed mesh generator. 
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INTRODUCTION 


Since the introduction of finite element analysis in the 
1960's and thru the 1 97 0 ' s it has become obvious that there was a 
clear need for powerful mesh generation techniques. 

What evolved were two different approaches to mesh gener- 
ation. One was to construct a model starting with building block 
elements and some computer aid to place and replicate these 
elements. The result was to build a model of the desired compo- 
nent or design. The second approach was to start with the compo- 
nent or design and break it first into topologically simple 
regions which in turn could be automatically broken into elements. 

In parallel with this program CAD systems were evolving. 

These systems provided powerful means for defining complex geometry 
and storing it in computer readable form. 

By 1978 we at Sikorsky were becoming more and more heavily 
involved in finite element modeling principally with manual 
methods. The burden of this was unacceptable. In our survey of 
available mesh generation systems we found no commercially avail- 
able system which we felt fully exploited the available tech- 
nology. 

It was our opinion that finite element modeling should 
naturally start from defined design geometry and follow the break- 
down method of mesh generation. Further we felt that the power 
of interactive graphics CAD systems should be fully exploited to 
obtain the geometry. 

We had available to us some elements of the desired system. 

We had available an in-house developed CAD system called CAIDS 
(Computer Aided Interactive Design System) and a mesh generator 
developed at Pratt & Whitney Commercial Products Division which 
met most of our requirements. With these elements we were well 
positioned to couple them into an effective system. 
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IFEMS CONFIGURATION & USE 


The total interactive modeling system is actually made up of 
three separate programs, (Figure 1). The first part is a complete 
computer aided design system, CAIDS, (Computer Aided Interactive 
Design System) modified to communicate with IFEMS. CAIDS repre- 
sents a highly sophisticated interactive design system equipped 
with a large software library of geometry generating routines. 

CAIDS hardware is the ADAGE (Figure 2) vector generator terminal 
with real time dynamic rotation, translation and zoom display 
features. The CAIDS display space is three dimensional, thus 
allowing for geometry definition of any conceivable component. 
Future features will include complete three dimensional math 
routines capable of solving for such things as intersections of 
general 3D surfaces. 

CAIDS uses 3D parametric rational cubic functions to store 
its geometry. These functions are stored as labeled groups or 
OBJECTS which are defined by the user. These OBJECTS become the 
means by which the analyst ultimately defines the finite element 
regions for the mesh generator and is the key to IFEMS. 

The second program is IFEMS. Its purpose is to accept the 
component geometry from CAIDS and combine it with the user de- 
fined model attributes (number of elements, boundary conditions, 
loads, etc.) and ultimately prepare an input file for the mesh 
generator. 

A functional overview of IFEMS is presented in Figure 3. 

The analyst begins by defining the actual component geometry using 
CAIDS. Upon completion of geometry definition, the analyst 
divides the model into regions. CAIDS then creates a neutral file 
which contains the model geometry grouped by OBJECTS. This neutral 
file is independent of the CAIDS data base and is used as a com- 
munications media to IFEMS. For large problems, the analyst may 
create several neutral files, the sum of which represent the 
compl ete model . 

Mesh generator regions are defined by using the OBJECT 
feature of CAIDS. That is, a generic name is defined by the user 
(region number and type, i.e. 25QUAD) to identify the region. 

The appropriate geometry defining that region is then assigned 
to the object. 
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Phase One of IFEMS is a preprocessor which reads the neutral 
file(s) performing four basic functions: (1) conducts basic data 

validation tasks on the region data looking for such errors as 
missing or non-adjacent edges, (2) reformats neutral file data 
into in-core arrays suitable for the downstream mesh generator 
(See Figure 4), (3) creates additional data, (from the cubic 
functions) required by the mesh generator such as node points 
representing the corners of the regions and (4) merges in previous- 
ly defined parts of the model. At the end of Phase One, a per- 
manent file, Perm One, is created (or updated in the case of 
multipart models) and stored on disk. The Perm One file contains 
all information generated during Phase One (See Figure 5). 

Phase One is an interactive processor (run on a Tektronix 4014 
terminal) and allows minor user intervention and minor corrections 
during data validation. 

Phase Two Part One of IFEMS is the major section. This is 
where the actual finite element model region attributes are 
defined: element type, number of elements per region, element 

spacing, material definition, boundary conditions, applied load, 
etc. This program begins by reading the region geometry from 
Phase One, via the Perm One File. 

The program uses a system of 40 menus to list the various 
options used in defining the region attributes (Figure 6). 

Through these menus the analyst controls the operational flow of 
the program. Displayed on the menu is the particular region that 
the analyst is working on. 

The menu system is hierarchical in that a given option 
automatically envokes a predetermined series of lower level menus 
(Figure 7 & 8). This system was designed to be tutorial. Given 
an upper level option, such as defining the number of elements 
along a region edge, the program will effectively guide the 
analyst by prompting him for the correct input. On the other hand 
the user can always interrupt the flow to return to the higher 
level menus. 

A Phase Two session would begin with the menu shown on 
Figure 9. The analyst will normally use Option D to define default 
region attributes, then Option E is used region by region to 
change only those attributes which are different from the default 
values. Changing the region attributes is accomplished by enter- 
ing the letter code corresponding to the desired change. Some of 
these changes are made to the current menu and others will auto- 
matically involve lower level menus (See Figures 7 & 8). 
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The user has additional control over the program through the 
use of the BAIL, CONTINUE, or NEXT options. BAIL nullifies all 
current input and reverts to the previously defined values. 
CONTINUE stores all current information and returns to the pre- 
vious menu for additional processing. NEXT is used to define 
the information for the next region or edge. 

Disk data sets are used to transfer information between pro- 
grams. The Perm Two file contains all the region attributes re- 
quired by the mesh generator. This file is updated each time the 
anlayst completes a region. After all the regions attributes 
have been defined the program will compress the Perm Two file and 
apply the default values to those regions which were not ex- 
plicitly accessed. 

Phase Two, Part Two is basically a pre-prpcessor to the mesh 
generator, PWAMESH . It creates an input file (from the Perm Two 
file) containing the region attributes which are then read by the 
mesh generator. PWAMESH then generates and displays the final 
mesh. The analyst reviews the final mesh and passes back and 
forth between Part 1 and 2 until he is satisfied. At this stage, 
PWAMESH is instructed to generate the bulk data required by the 
finite element program, typically NASTRAN. 

The third program, PWAMESH, is a general purpose three 
dimensional finite element mesh generator, developed by Pratt & 
Whitney Aircraft, Commercial Products Division. This program is 
executed from a Tektronix 4014 CRT. PWAMESH includes such 
features as variable thickness regions, generalized boundary 
condition definition, uncoupled regions for sliding mechanisms, 
automatic generation of quadratic temperature and pressure 
fields, automatic region joining as well as many other features. 

PWAMESH is also equipped with an interactive graphic package 
for reviewing the finite element model prior to running the finite 
element code. Its interactive review menu (Figure 10) will plot 
either the input coarse region(s) or the generated fine region(s) 
of the model. The entire model or particular regions can be dis- 
played at any viewing angle along with the element and node iden- 
tification numbers. Offline plots can be drawn using the VERSATEC 
plotting system. In addition, the user can blow-up any particular 
section by indicating the lower left and upper right boundaries 
of the model. The elements of the model can also be shrunk to 
make individual elements more visible (Figure 11). PWAMESH will 
prepare a complete finite element input deck for either NASTRAN, 
MARC or an in-house boundary integral program. 
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IFEMS SPECIAL FEATURES 


IFEMS was designed to minimize the amount of user defined 
input data. This has been achieved by allowing the user to define 
a set of initial model attributes. These attributes are applied 
to each region the first time the user displays the region. The 
user need only change those attributes that are unique for that 
region. Typically just the number of elements along an edge, or 
thickness. 

Additional features include automtic continuous backup of 
user input, free field input, automatic cursor positioning and a 
full compliment of displayable help pages. 


CONCLUDING REMARKS 


IFEMS, in conjunction with CAIDS and PWAMESH provide the 
analyst with a completely interactive tool for modeling mechanical 
component independent of geometric complexity. The usual errors 
associated with manually preparing mesh generator input are 
virtually eliminated as the analyst continually 'sees' the result 
of his work. 

Production use of IFEMS at Sikorsky has resulted in tremendous 
saving in both man-hours and lead time. Models that previously 
took three weeks are now being generated in two days. Solving 
complex problems and performing stress evaluation for minimum 
weight and cost with practical lead times has resulted in the more 
efficient use of the analyst's time and knowledge. 
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A TOTAL INTERACTIVE MODELING SYSTEM 



FIGURE 1 
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FIGURE 2 - ADAGE Vector Generator Terminal 
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PERM ONE FILE 
STRUCTURE 


REGION DEFINITIONS 

CORNER NODE DEFINITION 

COORDINATE SYSTEM DEFINITIONS 

REGION EDGE FUNCTION DEFINITION & POINTERS 

FUNCTION COEFFICIENTS 

ARRAY PARAMETERS (ie, NO. ENTRIES, ETC.) 

READ/WRITE UNIT NUMBERS 


FIGURE - 5 
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DEFINE ATTRIBUTES FOR REGION 1 (ENTER REGION ID THE FIRST TINE) 


SELECT ATTRIBUTE 

INPUT/CURRENT 


A-NO. ELS DIRECTION 1 / 6 

D-NO. ELS DIRECTION 2 / 6 

E-NO. ELS DIRECTION 3 / 6 

F-ELENENT TYPE CQUAD4 

G-ELENENT SPACING EQUAL 

J-EDGE INTERPOLATION 

((-MATERIAL ID / 1 

L-THICKNESS 1.000 


INPUT/CURRENT 

M-TEMPERATURE 0 . 0 

P-PRESSURE 0.0 

Q-FACE B.C. (TRANS) 

S-FACE B.C. (ROT) — 

T-SPECIAL FUNCTIONS 

U-IGNORE ELEMENTS 

U-REGION B.C. 

X-OUTPUT DISP SYS / 0 

Y-PERM CONSTAINTS / 0 


U-UIEU CHANGE 
N/C = NOT CONSTANT 


B-BAIL C-CONTINUE 
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ELEMENT SPACING FOR REGION 1 


QUAD REGION - PERMISSIBLE EDGES ARE 
1,2,3 OR 4 

ENTER EDGE NUMBER 


SELECT SPACING OPTION 


A UNEQUAL 
D GRADIENT 

E DISPLAY CURRENT UALUES 
F ENTER GRADIENT HERE 
G SAME AS EDGE OF REGION 
(NOTE: BE SURE THIS EDGE EXISTS, 
OTHERUISE PUAMESH UILL FAIL) 


SELECT 

B-BAIL 


U-UIEU CHANGE 



H-HELP N-NEXT EDGE 


C-CONTINUE 


FIGURE - 7 



ELEMENT TYPE FOR REGION 1 CURRENT ELEMENT TYPE: CQUAD4 

SELECT ELEMENT TYPE 

A CQUAD4 - 4 NODED ISOPARAMETERIC PLATE 

E CTRIA3 - 3 NODED UARIABLE THICKNESS TRIANGULAR PLATE, CONST STRAIN 

F CTRIAE - 3 NODED CONSTANT STRAIN, CONSTANT THICKNESS TRIA PLATE 

G CHEXA - 20 NODED SOLID, ANISOTROPIC/ISOTROPIC 

J CHEX20 - 20 NODED SOLID, ISOPARAMETRIC, ISOTROPC ONLY 

K CHEX8 - NOT AUAILABLE* USED ITEM 'l' 

L CHEXA - 8 NODED SOLID, ANISOTROPIC/ISOPTROPIC 

M CPENTA - 16 NODED ISOPARAMETERIC SOLID UEDGE 

P CPENTA - 6 NODED ISOPARAMETRIC SOLID UEDGE 

Q CBEAM - BEAM ELEMENT 

R CROD - AXIAL ROD 


B-BAIL C-CONTINUE H-HELP 
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CD 


SELECT ITEM 


ttt IFEPIS *** 

MAIN MENU 

OPTION PREVIOUSLY 
SELECTED 


D DEFINE REGION DEFAULT VALUES ( ) 

E DEFINE INDIVIDUAL REGION ATTRIBUTES C ) 

F DEFINE MATERIAL PROPERTY(S) < > 

G GENERATE PUAMESH INPUT < > 

J LIST REGION ATTRIBUTES FOR REGION - 
K GENERAL OPERATING INSTRUCTIONS 


T TERMINATE SESSION (UPDATE PERM 2 FILE) 

H HELP 

X DEBUG TRACE OFF 

U ABORT PHASE II 
V LIST COMPLETED REGIONS 
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SELECTS 

1 4-UERSATEC 

2 12-UERSATEC 

3 ADD NODE NO 

4 BLOU-UP 

5 BLOU-UP + 
NODE NO 

6 BLOU-UP + 
FIND NODE NO 

7 EXPLODED 
UIEU 

8 ADD ELM NO 

9 BLOU-UP + 
ELM NO 

10 OPTION 5+9 

11 CONTINUE 
ENTER- 


nOUBS - 10 


CO 

CO 



SELECT: 

1 ADD NODE NO 

2 BLOU-UP 

3 BLOU-UP + 
NODE NO 

4 BLOU-UP + 
FIND NODE NO 

5 ADD ELN NO 

6 BLOU-UP + 
ELM NO 

7 OPTION 3+6 

8 CONTINUE 
ENTER- 


X 0*0 Y 23.2 Z 34.3 
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APPLICATION OF A DATA BASE MANAGEMENT SYSTEM TO A 


FINITE ELEMENT MODEL 

James L. Rogers, Jr. 
Langley Research Center 


INTRODUCTION 


In today's software market, much effort is being expended on the develop- 
ment of data base management systems (DBMS). Most commercially available 
DBMS were designed for business use. However, the need for such systems with- 
in the engineering and scientific communities is becoming apparent. 

A potential DBMS application that appears attractive is the handling of 
data for finite element engineering models. The purpose of this paper is to 
explore the application of a commercially available, business-oriented DBMS 
to a structural engineering, finite element model. The model, DBMS, an 
approach to using the DBMS, advantages and disadvantages are explored in de- 
tail; and plans for research on a scientific and engineering DBMS are dis- 
cussed. 


THE FINITE ELEMENT MODEL 


Many organizations use more than one finite element computer program 
for analyzing structural models. For instance NASTRAN might be used because 
of its generality, standardization, and universal acceptance, or SPAR 
might be used to gain quick turnaround and interactive capability. Input 
into these programs usually consists of data representing a finite element 
model. Typically, the model has a single physical representation which 
can easily be pictured as a hierarchial structure composed of the complete 
model, substructures, elements, and grid points (see fig. 1). This single 
representation, however, can take the form of input to any one or more of 
the many finite element programs available today. 

If an engineer wishes to use more than one application program for 
analysis, a separate input deck must be maintained for each program, although 
in truth only one model is being analyzed. This proliferation of input 
decks can lead to a loss of integrity in the data, because of the difficulty 
in keeping all related input decks at the same level of modification — 
especially when many changes are being made to the model. To solve the 
proliferation problem and maintain integrity of the data, the finite element 
model input data can be stored in a data base. 
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THE DATA BASE MANAGEMENT SYSTEM 1 


Storing the model input data in a data base independent of all the 
application programs that are to make use of it as input can solve the 
problem of proliferation of input decks and aid in maintaining the integrity 
of the model. The data base is loaded into the DBMS in the steps shown 
in figure 2. Step 1 simply defines a data base name for the user to work 
with. The model input data must now be stored in data elements, the basic 
components of a data base. By using the data base definition the hierarchial 
structure of the finite element model is easily adapted to the hierarchial 
structure of the data base. A typical data base definition is shown in 
figure 3. Particular attention should be paid to its generality and its 
hierarchial arrangement. Most any finite element model, large or small, 
can be input with this definition. The hierarchial structure can be seen 
by scanning the indentation of the definition. Data element 1 allows 
input of a model name. Data element 2 is a repeating group (RG) of sub- 
structures. A repeating group allows more than one substructure per model 
to be input. Data element 3 allows input of a substructure name. Data 
element 4 is a repeating group of finite elements within a substructure. 

Data elements 5-10 provide for input of element information. Data element 

11 is a repeating group of grid points within a finite element. Data elements 

12-16 provide for input of grid point information. Data elements 20-49 

define supporting information (not necessarily hierarchial ) about the 

model such as element properties, single point constraints, omitted coordinates, 

and eigenvalue information. Step 3 (fig. 2) inputs the finite element 

model into the data base through the data base definition. 

After the data has been loaded into the data base, it cannot be used in 
an application program without undergoing some sort of conversion. The 
easiest approach when working with already existing programs, is to use 
FORTRAN pre-processor programs to convert the data from the independent 
data base format to an input format for a particular application program 
such as NAS TRAN or SPAR (see fig. 4). Pre-processor programs for converting 
data base input into NASTRAN and SPAR input decks are shown in Appendices 
I and II, respectively. A different pre-processor is needed for each 
application program that will use the model data as input. A data manipulation 
language (DML) which can be embedded within the pre-processor programs allows 
the user to write a FORTRAN program that will interact with the data base. 

The DML statements of the programs in the appendices have *PL in columns 1-3. 


SYSTEM 2000, a commercially available data base management system developed 
by MRI Systems Corporation, was used for all research done for this paper. 
However, almost any DBMS with a FORTRAN procedural language and query language 
could have been used. Use of commercial products and names of manufacturers 
in this report does not constitute an official endorsement of such products or 
manufacturers either expressed or implied, by the National Aeronautics and 
Space Administration. 
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Another powerful tool of some DBMS is the query language. This language 
allows the user to interactively query the data base. Using this language, 
he can list the grid points associated with a particular finite element or 
vice-versa he can list the finite elements associated with a particular grid 
point. Many other associations can be listed depending upon the relationships 
defined in the data base definition. 

This has only been a brief introduction into the world of DBMS. More 
details concerning the data base definition, data base input, DML, query 
language, and creation of a data base, can be obtained from manuals available 
from the developer of a specific DBMS. 


APPLILCATION OF A DBMS TO A FINITE ELEMENT CANTILEVER BEAM 


The data base input of a cantilever beam is shown in figure 5. Close 
attention should be paid to the correlation of numbers between figure 3 (data 
base definition) and figure 5 (data base input). To demonstrate the process 
shown in figure 4, a cantilever beam finite element model is created (see 
fig. 6). The data base definition and data base input described above and 
shown in figures 3 and 5 are used in this demonstration. After the data is 
stored in the data base, it is input into the pre-processor programs shown 
in appendices. Other input that does not pertain to the model but is necessary 
input to NASTRAN and SPAR (e.g., Executive and Case Control decks of NASTRAN) 
are read in from another file and inserted into the proper place in the output 
from the pre-processor programs. The pre-processor program for NASTRAN is not 
set up to handle PARAM and CNGRNT Bulk Data cards, thus these cards are also 
read from the same file as the Executive and Case Control cards. The pre- 
processor programs could be expanded and made sufficiently general to handle 
almost any input requirement. The pre-processor programs output a file ready 
for input into either NASTRAN (fig. 7) or SPAR (fig. 8). Thus it is shown 
that a finite element model data can be stored in a data base, independent of 
any application program which would use the model data as input. With the aid 
of pre-processor FORTRAN programs, the model data can subsequently be converted 
to a format for direct input into one of several application programs. 


ADVANTAGES AND DISADVANTAGES 


The advantages in using a data base to store a finite element model have 
been discussed in detail in previous sections so only a summary appears in 
this section. The data base requires only one copy of the model to be stored 
even though the data may be used as input to more than one application program. 
This situation aids in maintaining the integrity of the model , especially 
a model undergoing frequent changes. The stored data are independent of the 
application programs that use it, thus changes can be made to the programs 
and/or the data without having cross effects. A query language in the DBMS 
gives the user easy interactive access to the data, allowing listings of many 
types of data elements and cross-references between data elements. 


225 



The use of a DBMS is not without disadvantages. Of primary importance 
is the computer storage overhead (disk space) required to store a model. The 
small cantilever beam problem described earlier took almost n,000ig words of 
storage. Changing the NON-KEY values to KEY in the data base definition re- 
sults in the data base requiring 23,500-ig words of storage. (A KEY value per- 
mits the user to make sorts on that particular data element.) The user must 
judiciously choose which data elements are KEY to minimize storage overhead. 
Another disadvantage is that most business-oriented DBMS do not permit "E" 
formats (scientific notation). Since many structural terms involve very large 
or very small numbers, the "E" format is essential. One final disadvantage is 
the lack of support for matrix or vector input to the data base. Although 
this capability is not required to store finit. element model data, it would 
be particularly useful when storing and querying intermediate or final results 
from the application programs. 

These are the major advantages and disadvantages revealed during the course 
of the research. The user will have to make trade-offs in determining the most 
optimal way to use a DBMS with finite element models and computer programs. 


PLANS FOR DBMS WORK WITH FINITE ELEMENT MODELS 


New data base management systems (DBMS) are now being designed and de- 
veloped for engineering and scientific applications. For this reason, as well 
as the previously discussed disadvantages, no future research is planned 
using a commercial business-oriented DBMS with a finite element model. Future 
research is planned using E/S DMS, a scientific and engineering oriented DBMS 
being developed under a NASA grant at the University of Texas. E/S DMS is 
scheduled to be delivered to NASA Langley in August 1979. It executes on 
CDC CYBER computers with a NOS operating system. It has been proposed to imple 
ment E/S DMS on a minicomputer allowing users to distribute the data base as 
well as the processing between a mainframe and a minicomputer. This arrange- 
ment would give the user a very powerful tool; because the design and modeling 
could then be done on the faster (line speed) interactive graphics terminal, 
while the "number crunching" applications programs such as NASTRAN and SPAR 
could be executed on the more powerful mainframe. Future research is planned 
along these lines. 


CONCLUDING REMARKS 


A structural engineering, finite element model can be represented as a 
hierarchial structure within a business-oriented data base management system 
(DBMS). Use of the DBMS and conversion pre-processor computer programs, allow 
the data to be stored independent of application programs that use the data 
as input. This arrangement has the advantage that only one copy of the data 
needs to be stored and maintained even though more than one application program 
can use the data as input. The DBMS may also allow the user to query the data 
base as an aid in his modeling. The data manipulation language of the DBMS 
gives the application programmer the power and flexibility of FORTRAN combined 
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with the data storage capabilities of the data base. A small example problem 
demonstrates the use of a data base with one specific finite element model. 
Output of the conversion programs is shown in the form of input decks to 
both the NASTRAN and SPAR structural analysis application programs. Advantages 
and disadvantages are described to show that while although a DBMS created for 
business-oriented data can be useful for handling a finite element engineering 
model, it does not have all the capabilities required for engineering and 
scientific computing. 
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Figure 1.- Hierarchial structure of a finite element model 


STEP 1 - CREATE THE DATA BASE 


STEP 2 - LOAD THE DATA BASE DEFINITION 


STEP 3 - LOAD THE DATA BASE INPUT 



Figure 2.- Steps for loading a data base. 
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1* MODEL (NAME X (21 ) ) 5 
2* SUBSTRUCTURES (RG) ; 

3* SUBSTRUCTURE NAME (NAME X ( 1 0 ) IN 2)5 
4* SUBSTRUCTURE ELEMEnTS(RG IN 2)? 

5* ELEMENT ID (INTEGER NUMBER 9(7) IN 4) } 

6* ELEMENT NAME (NON-KEY NAME X(7) IN 4); 

7* XI VECTOR COMPONENT (NON-KEY DECIMAL NUMBER IN 4)5 

8* X2 VECTOR COMPONENT (NON-KEY DECIMAL NUMBER IN 4); 

9* X3 VECTOR COMPONENT (NON-KEY DECIMAL NUMBER IN 4); 

10* ELEMENT PROPERTY ID (INTEGER NUMBER 9(7) IN 4)? 

11* ELEMENT GRID POINTS (RG IN 4) ; 

12* GRID POINT ID (INTEGER NUMBER 9(7) IN 11)5 

13* X COORDINATE (NON-KEY DECIMAL NUMBER IN 11)5 

14* Y COORDINATE (NON-KEY DECIMAL NUMBER IN 11)5 

15* Z COORDINATE (NON-KEY DECIMAL NUMBER IN 11)5 

16* PERMANENT SPC (NON-KEY INTEGER NUMBER 9(7) IN 11)5 
20* ELEMENT PROPERTIES (RG IN 2)5 

21* PROPERTY ID (INTEGER NUMBER 9(7) IN 20)5 

22* MATERIAL ID (NON-KEY INTEGER NUMBER 9(7) IN 20)5 

23* AREA (NON-KEY NAME X ( 7 ) IN 20)5 
24* MOMENTS OF INERTIA (NON-KEY NAME X ( 7 ) IN 20)5 
25* YOUNGS MODULUS (NON-KEY NAME X(7) IN 20)5 
26* SHEAR MODULUS (NON-KEY DECIMAL NUMBER IN 20)5 
27* POISSONS RATIO (NON-KEY DECIMAL NUMBER IN 20)5 
28* MASS DENSITY (NON-KEY DECIMAL NUMBER IN 20)5 
29* PROPERTY NAME (NON-KEY NAME X(7) IN 20)5 
35* SINGLE POINT CONSTRAINTS ( RG IN 2)5 

36* SPC SET ID (NON-KEY INTEGER NUMBER 9(7) IN 35)5 
37* SPC GRID ID (NON-KEY INTEGER NUMBER 9(7) IN 35)5 
38* SPC COMPONENT NUMBER (NON-KEY INTEGER NUMBER 9(7) IN 35)5 
39* OMITTED COORDINATES (RG IN 2)5 

40* OMIT GRID ID (NON-KEY INTEGER NUMBER 9(7) IN 39)5 
41* OMIT COMPONENT NUMBER (NON-KEY INTEGER NUMBER 9(7) IN 39)5 
42* EIGENVALUE INFORMATION (RG IN 2)5 

43* EIGENVALUE SET ID (NON-KEY INTEGER NUMBER 9(7) IN 42)5 

44* METHOD (NON-KEY NAME X<7) IN 42)5 

45* LO FREQ (NON-KEY DECIMAL NUMBER IN 42)5 

46* HI FREQ (NON-KEY DECIMAL NUMBER IN 42) 5 

47* ESTIMATE OF ROOTS (NON-KEY INTEGER NUMBER 9(7) IN 42)5 

48* DESIRED ROOTS (NON-KEY INTEGER NUMBER 9(7) IN 42)5 

49* NORMALIZING METHOD (NON-KEY NAME X(7) IN 42)5 


Figure 3.- Data base definition for a finite element model. 
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Figure 4.- Flowchart of finite element data base process. 
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24®8. 0 - 
26® 7 ,0+10® 
27®. 3® 
28«86» 
29»PBAR® 
36«3b*20« 
37*1* 
38*26* 


39»40»13* 

41*6* 

39*40*2* 

4 1 ® 6® 
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Figure 5.- Data base input for a cantilever beam. 
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Figure 7.- NASTRAN input deck output from preprocessor program NTRNS2K. 
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gure 8. SPAR input deck output from pre-processor program SPARS 2K 
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APPENDIX A 


PRE-PROCESSOR PROGRAM TO CONVERT DATA IN DATA BASE TO NASTRAN INPUT FORMAT 


PROGRAM NTRNS2K(INPUT,OUTPUT*TAPE10=OUTPUT »TAPE9) 

C 

C THE PURPOSE OF THE PROGRAM IS TO CONVERT DATA FROM A DATA BASE 
C INTO NASTRAN BULK DATA FORMAT 

C 

INTEGER GRID(30),PSPC(30)»PID(10) , C5 * Cl 0 * C 12 , C21 * C22 * C 1 6 
INTEGER C36,C37,CaO*CA1,CA3,CA/,CA8 

INTEGER SChNME »RCODE *FILLER*PASSrtD,DATASN*DATASP»TlME,DATE»CYCLE* 

1 SEPSYM,ENTRYT, STATUS 
C 

C COMMON BLOCK FOR PICASSO DATA BASE 
C 

*PL C0MMBL0CK/PICASS0/SCHNME*RC0DE»FILLER*LASTDS*PASSWD»DATASN,DATASP, 
*PL LEVEL ♦ TIME ,DaTE i CYCLE . SEPSYM * ENTRYT » STATUS. 

C 

C SCHEMA FOR TOTAL STRUCTURE 
C 

*PL SCHEMA/LMONE OF PICASSO/C1. 

c 

C SCHEMA FOR REPEATING GROUP 2 (SUBSTRUCTURE NAMES) 

C 

*PL SCHEMA/LZERO OF PICASS0/C3. 

C 

C SCHEMA FOR REPEATING GROUP A (ELEMENTS) 

C 

*PL SCHEMA/LONE OF P IC ASS0/CB,C6» CT *CB*C9 ,C 1 0 . 

C 

C SCHEMA FOR REPEATING GROUP 11 (GRID POINTS) 

C 

*PL SCHEMA/LTWO OF P I CASSO/C 12 * C13 , Cl A , C 1 b , C 1 6 . 

C 

C SCHEMA FOR REPEATING GROUP 20 (ELEMENT PROPERTIES) 

C 

*PL SCHEMA/LTHRE OF PICASSO/C21 »C22 »C23»C2A * CEB » C26 * C2 7 * C20 » C2V. 

C 

C SCHEMA FOR REPEATING GROUP 36 (SINGLE POINT CONSTRAINTS) 

C 

*PL SCHEMA/LFOR OF PICASS0/C36»C37»C38. 

c 

c SCHEMA FOR REPEATING GROUP 39 (OMITTED COORDINATES) 

C 

*PL SCHEMA/LFIVE OF P ICASSO/CAO » CA 1 . 

C 

C SCHEMA FOR REPEATING GROUP A 2 (EIGENVALUE INFORMATION) 

C 

*PL SCHEMAZLATE OF PICASS0/CA3»CAA.CA5»CA6,CA7,CABtCA9. 

*PL END SCHEMAS. 

DIMENSION DECK I N ( 8 ) . XCORD(30) . YCORD(30) * ZCORD(30) 

DATA BEGBULK/10HBEGIN BULK/ 

C 

C FIRST READ IN THE EXECUTIVE AND CASE CONTROL DECKS AND STORE ON 
C T APE 1 0 

1 READ (9,2) (DECKIN(I) 1 1 = 1*8) 

2 FORMAT ( 8 A 1 0 ) 

KRITE(10»2) ( DECK I N ( I ) *1 = 1,8) 

IFIDECK IN ( 1 ) ,NE . BEGBULK ) GO TO 1 
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c 

*Pl START S2K. 

PASSWD=6HR0GERS 

c 

C OPEN PICASSO DATA BASE 
C 

*PL OPEN PICASSO. 

C 

C CHECK RETURN CODE 
C 

IF (RCODE.EG. 0) 60 TO 10 
PRINT S.RCODE 

5 FORMAT ( # DATA BASE CAN NOT BE OPENED - RETURN CODE = 
60 TO 9 00 
C 

C CHECK STATUS 
C 

10 IF (STATUS. EQ.O) 60 TO 20 
PRINT IS. STATUS 

lb FORMAT (o DATA BASE DAMA6ED - STATUS = *,I2) 

60 TO 900 
C 
C 

C CREATE ELEMENT CARDS 
C 

20 RE AD ( 9 . 25 ) C3 
25 FORMAT(AIO) 

IF (EOF ( 9 ) ) 800,30 
30 CONTINUE 
1 = 0 
M = 0 
C 

C RETRIEVE SUBSTRUCTURE FROM DATA BASE 

C 

C 

*Pl 6ET1 LZERO WHERE C3. 

C 

C CHECK RETURN CODE 
C 

IF (RCODE.EQ. A) 60 TO 700 
IF (RCODE.EQ.O) 60 TO 39 
PRINT 35.RC0DE.C3 

35 FORMAT (<* RETURN CODE = *»I3.« FOR SUBSTRUCTURE ».A10> 
60 TO 900 
39 ISW=0 
C 

AO CONTINUE 
ISW=ISW*1 
C 

C RETRIEVE ELEMENT IDS FOR SUBSTRUCTURE FROM DATA BASE 
C 

*P L 6ETD LONE NEXT. 


13) 
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C CHECK RETURN CODE 
C 

IF (RC0DE.EQ.4) GO TO 49 
IF (RCOOE.Ey.O) GO TO 50 
PRINT 45 « RCOOE * C3 

45 FORMAT (» RETURN CODE = <**I3»* FOR SUBSTRUCTURE **A10** AT CHKPT2«) 
GO TO 900 

49 ISW=-l 

50 CONTINUE 
C 

C GET ELEMENT GRID INFORMATION 
C 

*PL GETO LTWO NEXT. 

C 

C CHECK RETURN CODE 
C 

IF (RC0DE.EQ.4) GO TO 70 
IF(RCODE.EQ.O) GO TO 60 
PRINT 55 * RCODE * C5 

55 FORMAT (* RETURN CODE = **I3.» FOR ELEMENT «*IT) 

GO TO 900 
60 1=1+1 

GRID(I)=C12 
XCORD < I ) =C 1 3 
YCORD ( I ) =C 1 4 
ZCORD ( I ) =C16 
PSPC<I)=Clb 
IF(ISW.EQ.l) GO TO 81 
C 
C 

C WRITE ELEMENT INFO ON UNIT 10 
C 

70 WRITE(10»80) CM6*M5»M10»GRID ( 1-1 ) * GR I D < I ) *CM/*CM8iCM9 

80 FORMAT (A8*4I8»3F8.1»7X* <t l <t *8X) 

81 M=M* 1 
CM6=C6 
M5=C5 
M10=C10 
CM7=C7 
CM8=C8 
CM9=M9 
PID (M) =C 1 0 

IF ( ISW.NE.-l ) GO TO 40 
C 

C REORDER GRID POINTS IN INCREASING OROER 
C 

85 DO 9b K=2*I 
LP=I-K*1 
DO 90 L=1*LP 

IF(GRIO<L>.LE.GRID(L+l> ) GO TO 90 
1 TEMP = GR I U ( L * 1 ) 

GHID(L*1)=GRID(L) 

GRID (L ) = I TEMP 
TEMP=XCORD(L*l) 

XCORO(L*l)=XCORD(L) 

XCORD(L) =temp 
TEMP=YC0RD(L+1) 

YCORD (L* 1 > = YCORD (L ) 

YCORO(L> = TEMP 
TEMP=ZCORD(L*l) 

ZCOHD(L*I > =ZCORD(L> 

ZCORO(L) = temp 

1TEMP=PSPC < L ♦ 1 > 

PSPC(L+1)=PSPC(L) 

PSPC (L) =ITEMP 
90 CONTINUE 
95 CONTINUE 
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c 

C WRITE GRID DATE ON UNIT 10 
C 

ITEmP=0 

uo no j=i«i 

IF(ITEMP.EQ.GRIO(J) ) GO TO UO 

WRITE ( 10, 100) GRID(J) *XCORD(J> ,PSPC(J) 

100 F ORM AT ( *GR 1 D * , 1 8 , BX » F8 . 2 , 24X , 1 8, 1 6X ) 

I TEMP=GR I 0 ( J ) 

110 CONTINUE 

c 

C REORDER PROPERTY IDS INTO INCREASING ORDER 
C 

DO 114 I = 2,m 
UP = M-I+l 
DO 111 L= 1 , LP 

IF(PID(L).LE.PID(L*1)) GO TO 111 
1PI0 = P I D (L ♦ 1 > 

PI D (L* 1 ) = PID(L) 

P ID (L > = IPID 
111 CONTINUE 
114 CONTINUE 
C 

*PL 6ET1 LTHRE FIRST. 

C 

C RETRIEVE ELEMENT PROPERTIES 
C 

IPID - 0 
C 

DO 150 1=1, M 

IF (PID(I) .EQ.IPID) GO TO 150 
C 

C21 = PID(I) 

4PL GET1 LTHRE WHERE C21. 

C 

C CHECK RETURN CODE 
C 

IF (RCODE.EQ.O) GO TO 120 
PRINT lib , RCODE , P ID ( I ) 

115 FORMAT (» RETURN CODE = *,I3,« FOR PROPERTY ID *,I7> 
GO TO 900 
C 

C WRITE MATERIAL PROPERTIES ON UNIT 10 
C 

120 WRITE(10*130) C29,C21,C22,C23,C24 

130 FORMAT (AG, 218, 2AG,40X) 

WRITEU0.140) C22,C2b,C27,C2B 
140 F ORM AT ( # M AT 1 * , I B , AB , BX , 2F8 . 2 , 32X ) 

IPID = PIO (I) 
lbO CONTINUE 
*PL GET! LFIVE FIRST. 

GO TO lbl 

160 CONTINUE 
C 

C RETRIEVE OMITTEO GRID POINTS 
C 

*PL geti lfive next. 

c 

C CHECK RETURN CODE 
C 

161 IF(RCODE.EO.O) GO TO 170 
IF (RCODE. EQ. 4) GO TO 190 
PRINT 16b, RCODE 

165 FORMAT (* RETURN COOE FOR LFIVE IS ««I3) 

GO TO 900 
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c 

C WHITE OMITTED COORDINATES ON Unit 10. 

c 

1T0 WHITE ( 10. 180) 0*0.041 
180 FORMA T ( *0M I T *,2I8,36X) 

go ro loo 

190 CONTINUE 
«PL 0ET1 UFOR FIRST. 

GO TO 102 

101 CONTINUE 
C 

C RETRIEVE SINGLE POINT CONSTRAINTS 
C 

*PL GET! LFOR NEXT. 

C 

C CHECK RETURN CODE 
C 

102 IF (RCODE. EG. 4) GO TO 220 
IF (RCODE.EJ. 0) GO TO 200 
PRINT 10s. RCODE 

195 FORMAT!* RETURN CODE FOR LFOR IS ».I3> 

GO TO 900 
C 

C WRITE SINGLE POINT CONSTRAINTS ON UNIT 10 
C 

200 WRITE ( 10,210) C36 , C37 , C38 
210 FORMAT t *SPC *«3IB»48X) 

GO TO 101 
220 CONTINUE 
»PL GET1 LATE FIRST. 

GO TO 224 
221 CONTINUE 
C 

C RETRIEVE EIGENVALUE INFORMATION 
C 

*PL GET 1 LATE NEXT. 

C 

C CHECK RETURN CODE 
C 

224 IF (RCODE. EO.O) GO TO 230 
PRINT 22S, RCODE 

22S FORMAT!* RETURN CODE FOR LATE = *.I3> 

GO TO 000 
C 

c write eigenvalue information on unit 10 
c 

230 WRI TE ( 10,240 ) C43,C44,C4S,C4b«C4 l , C4B 
240 F ORMAT ( *E I GR * . I 8 , A8 , 2F8 . 1 , 218 < 1 bX . * , BC » , 3X ) 

WRITE ( 10,230) C49 
250 FORMAT (*+dC «,Ad,64X) 

C 

C REAO REST OF OECK IN ANO STORE ON UNIT 10 
C 

252 READ(9.2) (DECKIn(I) ,1=1,8) 

IF (EOF (9) ) 255,253 

253 WRir£<10.2) (OF.CKIN(I) ,1=1,8) 

GO TO 252 

C 

C WRITE ENDDATA ON UNIT 10 
C 

255 WRITE ( 10, 2t>0) 

260 FO RMAT(*ENDDATA »,72X) 

ENUFILE 10 
REWIND 10 
GO TO 999 
700 PRINT 750. C3 

750 FORMAT!* NO SUBSTRUCTURE WITH NAME *,A10> 

GO TO 20 
BOO PRINT 850 

830 FORMAT (» END OF DATA INPUT*) 

C 

C WRITE FINAL ERROR MESSAGE 
C 

900 PRINT 950 

930 FORMAT!//* PROGRAM TERMINATED FOR ABOVE REASON*//) 
999 STOP 

End 
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APPENDIX B 


PRE-PROCESSOR PROGRAM TO CONVERT DATA IN DATA BASE TO SPAR INPUT FORMAT 


PROGRAM SPARS2K ( INPOT , OUTPUT * TAPER v TAPE 1 O-OUTPUT ) 

DIMENSION DECK IN (8) tlSPC(fa) ,IGRD(4) 

INTEGER CS»Ci0»C12tCEl*C22»Cib»C3b*C3(»C40*C4l»C43.C47tC48*C18 
INTEGER SCHNME tRCODE .FILLER .PASSED* DAT ASN . DAT ASP « T I ME * DATE . CYCLE. 

1 SEPSYMtENTRYT»STATUS 

c 

C COMMON BLOCK FOR PICASSO DATA BASE 
C 

*PL C0MM6L0CK/P I CASSO/SCHMNE.RCODE.FILLER.LASTDS .PASSwD.DAT ASN .DATASP. 

*PL LEVEL *TIME .DATA. CYCLE. SEPSYM. ENTRY!, STATUS. 

C 

C SCHEMA FOR TOTAL STRUCTURE 
C 

*PL SCHEMA/LMONE OF PICASSO/C1. 

c 

c SCHEMA FOR REPEATING GROUP 2 (SUBSTRUCTURE NAMES) 

C 

«PL SCHEMA/LZERO OF PICASSO/C3, 

C 

C SCHEMA FOR REPEATING GROUP 4 (ELEMENTS) 

C 

»PL SCHEMA/LONE OF P I CASS0/C5 . Cb » C 7 » C8 . C9 . C 1 0 . 

C 

C SCHEMA FOR REPEATING GROUP 11 (GRID POINTS) 

C 

«PL SCHEMA/LTWO OF PICASS0/C12»C13.C14»C13,C16„ 

c 

C SCHEMA FOR RtPEATING GROUP 20 ( ELEMENT PROPERTIES) 

C 

*PL SCHEMA/LTHRE OF PI CASS0/C2 1 , C22 . C2 3 ♦ C24 . C2S * C2b » C2 ( . C28 , C29 . 

c 

C SCHEMA FOR REPEATING GROUP 3b (SINGLE POINT CONSTRAINTS) 

C 

*PL SCHEMA/LFOR OF PICASS0/C3b,C37 ,C38. 

C 

C SCHEMA FOR REPEATING GROUP 39 (OMITTEO COORDINATES) 

C 

#PL SCHEMA/LFIVE OF PICASSO/C40,C41. 

c 

c SCHEMA FOR REPEATING GROUP 42 (EIGENVALUE INFORMATION) 

C 

#PL SCHEMA/LATE OF PICASS0/C43»C44«C4S.C4b,C47»C4tt,C49. 

*PL END SCHEMAS. 

C 

DATA CONMAT/1 OHMATER I AL C/ 

C 

C READ IN DESCRIPTIVE INFORMATION ANO CARDS THROUGH MATERIAL CONSTANTS 
C 

1 READ(9»2) ( DECK I N ( I ) . I = 1 . 8 ) 

2 FORMAT (8A10 I 

WRITE(10,2) (DECKIN(I) »I=1*8) 

1F(DECKIN(1) .NE.CONMAT) GO TO 1 
*PL START S2K. 

PASSWD=bHHOGERS 

C 

C OPEN PICASSO DATA BASE 
C 

*PL OPEN PICASSO. 
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C CHECK RETURN CODE 
C 

IF(RCODE.EQ.O) GO TO 10 
PRINT 5 . RCODE 

5 FORMAT (* DATA BASE CANNOT BE OPENED - RETURN CODE = *.13) 
GO TO 900 

CHECK STATUS 

10 IF (STATUS, EQ.O) GO TO 20 
PRINT 15.STATUS 

15 FORMAT ( * DATA BASE DAMAGED - STATUS = *.I2> 

GO TO 900 

RETRIEVE STRUCTURE FROM DATA BASE 

20 CONTINUE 
PL GET 1 LMONE NEXT. 

CHECK RETURN CODE 

IF (RCODE. EQ.O) GO TO AO 
PRINT 35 » RCODE 

35 FORMAT (* RETURN CODE = *.13.* FOR STRUCTURE* ) 

GO TO 900 
AO CONTINUE 

RETRIEVE MATERIAL PROPERTIES FROM DATA BASE 

PL GETD LTHRE NEXT. 

CHECK RETURN CODE 

IF (RCODE. EQ.O) GO TO 50 
PRINT A5. RCODE 

A 5 FORMAT (* RETURN CODE = *.13,* FOR MATERIAL CONSTANTS*) 

GO TO 900 

WRITE MATERIAL CONSTANTS ON UNIT 10 

50 WR I TE ( 1 0 * 55 ) C25.C27.C2B 
55 FORMAT ( * 1 * ♦ A7 » 1 X * F3. 2 , 1 X • F3 . 2 ) 

GET JOINT INFORMATION FROM DATA BASE 

WRITE<10»61) 

61 FORMAT (* JOINT LOCATIONS*) 

60 ISW=0 

62 CONTINUE 

PL GETD LONE NEXT. 

CHECK RETURN CODE 

IF (RCODE. EQ. A) GO TO 69 
IF (RCODE. EQ.O) GO TO 70 
PRINT 65. RCODE 

65 FORMAT ( * RETURN CODE FOR LONE = *.I3) 

GO TO 900 

69 ISW=-1 

70 CONTINUE 

*PL GETD LTWO NEXT, 
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CHECK RETURN CODE 

73 IF (RCODE. EQ. 4) GO TO 62 
IF < RCODE , EQ , 0 ) GO TO 60 
PRINT 7b. RCODE 

75 FORMAT < * RETURN CODE FOR LTtfO = *.I3 
GO TO 900 

WRITE JOINT LOCATIONS ON UNIT 10 

80 WRITE (10*65) C12.C13 
85 FORMAT (I2.1X.F3.1) 

IF < ISW.NE.-l ) GO TO 62 
90 CONTINUE 

RETRIEVE CONSTRAINTS FROM THE DATA BASE 
WR I TE ( 1 0 * 95 ) 

95 FORMAT ( * CONSTRAINT CASE 1*) 
ISPC(l)=Clb/10G0Q0 
DO 110 I = 2.6 
K=6-I 

ISPC(I)=C16/10«*K 
N= I - 1 
M= I 

DO 105 L = 1 ♦ N 
M-M- 1 

ISPC(I)=ISPC<I)-<ISPC(L>*10##M> 

105 CONTINUE 

110 CONTINUE 

DO 111 I = 1.6 
IF(ISPC(I) .NE.0) GO TO 112 

111 CONTINUE 

112 IF ( I . GE • 6 ) GO TO 120 
1 * 1-1 

K=6-I 

DO 115 J = 1 » K 
ISPC(J)*ISPC(J+I) 

115 CONTINUE 

120 WRITE (10.125) ( ISPC(J) . J=1 .K) 

125 FORMAT ( * ZERO #.6{I1»1X)) 

WRITE (10.130) 

130 FORMAT ( * * XQT ELD#) 

WRITE (10.135) 

135 FORMAT (#E23#/#NSECT=1#) 

HO CONTINUE 
PL GETD LONE FIRST. 

GO TO 1 A3 
142 CONTINUE 

GET ELEMENTS 

#PL GETD LONE NEXT. 



c 

C CHECK RETURN CODE 

c 

143 IF (RC0DE.EQ.4) 60 TO 200 
IF (RCODE.EO, 0) GO TO 150 
PRINT 145,RC0DE 

145 FORMAT ( # RETURN CODE FOR LONE = * * I 3 ) 
60 TO 900 
150 CONTINUE 
C 

C GET GRID POINTS FOR ELEMENTS 
C 

DO 170 1=1*4 
*PL GETD LTWO NEXT. 

C 

C CHECK RETURN CODE 
C 

IF (RC0DE.EQ.4) GO TO 180 
IF (RCODE.EQ.O) GO TO 160 
PRINT 156* RCODE 

155 FORMAT!# RETURN CODE FOR LTWO = #»I3) 
GO TO 900 
160 IGRD ( I ) =C12 
170 CONTINUE 
C 

C WRITE GRID POINTS ON UNIT 10 
C 

180 1=1-1 

WRITE (10, 190) (IGRD(J) *J=1»I) 

190 FORMAT (4(12, IX) ) 

GO TO 142 
C 

C WRITE EXIT 
C 

200 WRITE ( 10,210) 

210 FORMAT ( # • XQT EXIT#) 

GO TO 999 


C 

C ERROR EXIT 
C 

900 PRINT 910 

910 FORMAT (/# JOB TERMINATED DUE TO ABOVE ERROR#) 
999 CONTINUE 
C 

C CLOSE PICASSO 
C 

#PL CLOSE PICASSO. 

#PL END PROCEDURE. 

STOP 

END 
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